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The  purpose  of  this  thesis  was  to  investigate  a  method 
of  extending  the  range  of  identification  of  nuclear  sources. 
The  method  was  centered  on  Bayes*  theorem,  and  the  results 
were  compared  with  the  method  of  photopeak  identification, 
which  is  perhaps  the  most  common  method  used  today.  The 
Finite  Element  Method  was  used  to  transport  the  source 
spectra  out  to  various  ranges  and  Bayes*  theorem,  in 
conjunction  with  various  statistical  distributions,  was  used 
to  analyze  the  results. 

I  would  like  to  thank  Dr.  Larry  McKee  for  sponsoring 
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and  the  Harris  800  operating  system.  Without  him,  none  of 
this  work  would  have  been  possible.  In  addition,  I  must 
express  my  gratitude  to  Mr.  Seth  Tuuri,  for  asking  the 
fundamental  questions  that  keep  basic  research  alive. 
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Abstract 


This  research  effort  studies  the  application  of  Bayes' 
decision  theorem  to  extending  the  range  of  remote 
identification  of  nuclear  materials.  The  Finite  Element 
Method  was  used  to  develop  a  radiation  transport  code  which 
was  used  to  reconstruct  a  group  of  three  sample  sources  at 
distances  of  50,  100,  200,  300,  400  and  500  meters.  Both  the 
Poisson  and  multinomial  distributions  were  then  used  to 
simulate  measured  sources  in  a  low  count  environment  at  these 
six  ranges.  Bayes'  theorem  was  applied  to  the  resulting 
measured  sources  to  test  for  positive  identification. 

The  results  show  that  a  low  resolution  detector  can 
increase  the  range  of  remote  detection  an  average  of  100 
meters  when  compared  with  the  method  of  photopeak 
identification.  Bayes'  theorem  is  unable,  however,  to 
identify  sources  not  contained  in  the  library  of  known 


sources . 
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IDENTIFICATION  OF  NUCLEAR  MATERIALS  FROM  REMOTE  DETECTION  OF 

CHARACTERISTIC  GAMMA  RAYS 


I .  Introduction 


Background 

This  thesis  addresses  the  problem  of  identifying 
specific  radioactive  sources  by  examining  their 
characteristic  gamma  ray  spectra.  The  primary  effort  will  be 
to  develop  a  technique  of  extending  the  range  of  reliable 
detection  of  nuclear  materials. 

Normally,  the  analysis  of  characteristic  gamma  ray 
spectra  involves  matching  the  photopeaks  of  the  given 
spectrum  with  those  present  in  the  spectra  of  known  isotopes 
for  the  purpose  of  identifying  the  constituents  of  the 
source.  Such  identification  of  nuclear  material  is  of 
interest  in  order  to  learn  the  nature  of  given  sources,  that 
is,  to  identify  the  'ingredients'  of  a  nuclear  source. 
However,  current  methods  which  use  photopeak  analysis  for 
identification  purposes  are  of  limited  value  once 
downscattering  of  the  peak  energies  has  occurred. 

In  order  to  be  of  practical  benefit,  the  identification 
of  nuclear  materials  must  consider  the  non-ideal  case  of 
remote  detection.  In  other  words,  the  given  gamma  ray 
source  spectrum  must  be  examined  at  a  distance  and  therefore, 
will  have  been  transported  through  the  atmosphere  for  a 
predetermined  distance.  Transporting  such  a  spectrum  through 


air  introduces  a  number  of  complications  due  to  the 
interaction  of  photons  with  molecules  and  atoms  in  the 
earth's  atmosphere. 

By  interacting  with  the  particles  of  the  air,  photon 
energy  is  degraded  by  three  main  processes:  the 
photoelectric  effect,  compton  scattering,  and  pair 
production.  The  importance  of  each  process  depends  on  the 
initial  energy  of  the  photon.  It  is  important  here  to 
realize  that  the  compton  scattering  process  is  responsible 
for  degrading  the  ideal  (untransported)  gamma  spectrum  by 
'smearing'  the  energy  of  the  photons  over  a  wide  range  of 
values.  Therefore,  Compton  scattering  is  the  primary 
mechanism  for  the  degradation  of  gamma  photopeaks  in  this 
problem . 

The  Problem 

This  thesis  will  study  the  ability  to  discriminate 
between  different  gamma  ray  source  spectra  by  using  Bayes' 
theorem.  Bayes'  theorem  will  be  employed  as  a  decision  aid 
to  assess  the  probability  of  having  a  certain  source,  given 
the  transported  spectrum.  By  plotting  the  probabilities  as  a 
function  of  distance  from  the  source,  we  can  determine  the 
reliability  of  Bayes'  theorem  as  a  correlation  technique  for 
extending  spectral  identification  beyond  the  current 
limitations  of  photopeak  analysis. 

Simply  stated,  this  thesis  will  test  the  assumption  that 
Bayes'  theorem  will  work  at  distances  from  the  source  where 
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less  than  ideal  circumstances  include  measuring  a  given 
source  once  at  extreme  range  with  a  small  number  of  counts 
collected  . 

Methodology:  the  Finite  Elements  Method 

The  method  used  to  approximate  the  solution  to  the 
diffusion  equation  is  the  Finite  Elements  Method,  It  was 
chosen  because  of  its  accuracy  in  comparison  to  other  methods 
such  as  the  Method  of  Integration  and  it  is  at  the  same  time 
less  complicated  to  program  than  more  advanced  methods  (8). 

In  general,  the  diffusion  equation  is  expressed  as 

-  7«(DVF)  +  IaF  -  S  (1) 

where  the  variable  D  is  the  diffusion  coefficient,  E  is  the 
macroscopic  absorption  cross  section,  S  is  the  source 
function,  and  F  is  the  fluence.  Also,  it  is  important  to 
note  that 
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where,  Etis  the  macrosopic  transport  cross  section.  The 
Finite  Elements  Method  assigns  a  penalty  function  given  by 


P(F) 


i[dV(DVF»VF  +  ZaF  2 
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The  radiation  transport  problem  will  be  solved  by 
implementing  a  numerical  solution  to  the  diffusion  equation. 
The  diffusion  equation  will  then  be  used  to  solve  for  the 
scattered  fluence  contribution  of  each  energy  group  in  each 
source.  This  scattered,  or  diffuse  solution  (including 
downscattering  and  inscattering  contributions)  will  be  added 
to  the  uncollided  fluence  (virgin  fluence)  solution  of  each 
energy  group,  resulting  in  a  value  of  total  fluence  for  each 
energy  group. 

By  solving  the  diffusion  equation  for  each  library 
source  at  a  number  of  different  ranges,  the  reference  sources 
are  arrived  at,  which  will  be  discussed  later  in  this 
chapter.  As  the  distance  increases  from  the  source,  the 
uncollided  fluence  drops  off  exponentially  (see  appendix  A), 
and  at  some  point,  there  is  an  extremely  small  contribution 
to  the  solution  due  to  the  uncollided  fluence,  but  still  a 
sizable  contribution  from  the  diffuse  (collided)  fluence.  In 
other  words,  the  further  away  from  the  source  the  measurement 
is  taken,  the  larger  is  the  percentage  of  diffuse  fluence. 

At  some  very  large  distance  then,  all  the  photons  detected 
will  have  been  scattered  at  least  once,  and  therefore  no  more 
uncollided  fluence  is  present.  The  concept  of  using  Bayes' 
theorem  as  an  analysis  tool  is  based  on  using  the  diffuse 
solution  in  addition  to  the  uncollided  solution  to  identify 
sources  measured  under  less  than  ideal  circumstances.  These 


In  the  last  chapter,  a  number  of  source  spectra  were 
introduced,  based  on  high  resolution  data  obtained  at  the 


surface  of  each  source.  Recall  that  the  spectral 
distribution  for  each  source  was  given  in  counts  per  gamma 
photon  per  unit  energy,  or  alternately,  in  counts  per  photon 
per  channel  (each  energy  group  representing  a  channel).  In 
this  chapter,  the  methodology  will  be  developed  to  transport 
each  of  these  library  spectra  to  any  given  distance  away  from 
the  original  source.  The  results  of  the  radiation  transport 
of  each  of  the  three  library  sources  will  be  presented  at  the 
end  of  this  chapter. 

Assumptions 

In  order  to  simplify  the  radiation  transport  problem 
into  a  manageable  one,  the  following  assumptions  are  made: 

(1)  The  source  photons  are  transported  through 
homogeneous  air. 

(2)  The  source  and  the  detector  are  located 
at  sea  level. 

(3)  An  ideal  detector  is  assumed,  so  that  all  photons 
reaching  the  detector  are  counted. 

(4)  The  fluence  of  photons  is  radially  symmetric 
(no  angular  dependence). 
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TABLE  I 1-2 


Relative  Intensity  Distribution  of  Library  Sources 

(counts  /  photon) 

Group  Number  Source  A  Source  B  Source  C 


2.20E-4 

2.37E-4 

1.47E-4 

7.92E-4 

9.06E-4 

8.35E-4 

1.70E-3 

5.46E-4 

1.43E-3 

2. 11E-3 

3.09E-4 

0.0 

6.99E-3 

3.24E-4 

2.27E-3 

1.00E-2 

7.07E-4 

0.0 

2.47E-2 

1.35E-2 

2.08E-3 

5. 16E-1 

7.90E-2 

1.28E-2 

3.06E-1 

1.91E-1 

1.28E-2 

4.4&E-2 

1.87E-1 

1 .83E-1 

5.09E-2 

3.24E-1 

4.43E-1 

2.58E-2 

1.34E-1 

2.41E-1 

9.80E-2 

6.73E-2 

1.00E-1 

TABLE  I I- 1 


"detector"  is  provided  in  table  II-l.  The  18  gamma-group 
cross  section  data  file  is  attached  in  appendix  B. 

The  second  step  involved  totalling  the  counts  contained 
in  each  of  the  newly  defined  eighteen  channels  for  each 
source.  In  addition,  the  background  data  was  subdivided  into 
the  eighteen  channel  structure. 

Finally,  the  background  was  subtracted  from  each  of  the 
sources,  and  the  total  number  of  counts  was  determined  for 
each.  The  relative  intensity  of  each  source  was  determined 
by  dividing  the  number  of  counts  in  each  channel  by  the  total 
number  of  counts  for  that  particular  source.  These  relative 
intensities  are  provided  in  Table  II-2.  The  final  library 
spectra  are  pictured  in  figure  1,  and  a  listing  of  the 
relative  intensities  in  counts  per  KeV  per  photon  are  given 
in  table  II-3. 

As  can  be  seen  in  figure  1,  all  three  sources  have 
identifying  features  that  makes  each  one  unique.  Sources  B 
and  C  show  the  most  similarities  with  respect  to  each  other, 
and  can  be  used  as  a  benchmark  during  the  analysis  section, 
as  mentioned  earlier.  As  the  transport  distance  increases, 
downscattering  of  gamma  photons  should  make  sources  B  and  C 
'converge'  to  a  common  spectrum,  whereas  source  A  should 
maintain  its  unique  identity  for  larger  ranges. 

In  the  next  chapter,  the  results  of  the  radiation 
transport  code  development  will  be  used  to  generate  the 
measured  and  reference  sources,  based  on  the  library  sources 
presented  in  this  chapter. 


photopeak  identification  requires  that  a  photopeak  be  at 
least  three  standard  deviations  above  the  background  level 
for  positive  identification  (three  sigma  being  99.7  percent 
probability  of  certainty).  Please  reference  appendix  A  for 
more  details  on  this  subject. 

Some  notable  photopeaks  observed  in  the  high  resolution 
specta  include  the  767  KeV  and  1001  KeV  lines  from  238U  and 
the  375  KeV  and  414  KeV  lines  of  239Pu.  Numerous  other 
photopeaks  from  other  radioisotopes  contribute  to  the  overall 
relative  intensities  in  each  channel. 

The  second  criterion  listed  above  will  test  how  Bayes 
Decision  Theory  can  discriminate  between  two  similar  nuclear 
sources,  in  contrast  to  two  sources  which  have  dissimilar 
spectral  features.  In  addition,  a  fourth  source  will  be 
introduced,  which  will  be  different  from  all  of  the  library 
spectra,  in  order  to  examine  the  effects  of  an  "unknown" 
source  on  the  analysis. 

Creation  of  the  Low  Resolution  Sources. 

Having  selected  the  three  high  resolution  spectra  (in 
line  w  th  the  above  criteria),  the  task  remains  to  translate 
these  into  the  three  low  resolution  library  spectra  for  use 
in  this  problem.  This  was  accomplished  in  three  basic  steps. 

First,  the  high  resolution  spectral  data  was  divided 
into  an  eighteen  energy  group  structure,  which  was  used  in 
thii  problem.  The  eighteen  group  structure  was  chosen 
because  the  gamma  photon  cross  section  data  was  readily 
available  (14).  The  energy  structure  of  this  low  resolution 
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II .  Source  Selection 


The  first  order  of  business  in  this  problem  is  to 
identify  several  source  spectra  that  will  represent  the  gamma 
ray  signature  of  different  nuclear  materials.  These  source 
spectra,  or  "library  spectra,"  will  then  be  used  to  test  the 
feasibility  of  Bayes  Theorem  as  a  remote  identification  tool. 
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The  Library  Spectra 

Method  of  Selection. 

There  are  two  possible  approaches  to  the  source 
selection  process.  In  the  first  approach,  the  library 
spectra  are  chosen  to  represent  the  gamma  spectra  of  specific 
radioisotopes,  that  is,  each  source  is  composed  of  a  simple 
combination  of  radioisotopes,  with  the  intensity  levels 
determined  arbitrarily.  Using  a  second  approach,  the  library 
spectra  represent  actual  nuclear  materials.  By  using  the 
second  method,  the  intensity  levels  of  the  library  spectra 
are  more  realistic,  having  been  taken  from  actual 
measurements  of  different  nuclear  materials. 

The  second  method  of  source  selection  was  used  in  this 
research  project,  in  order  to  more  realistically  assess  the 
utility  of  Bayes  Theorem  as  a  method  of  remote  source 
identification. 

Spectral  Data. 

A  number  of  measured  gamma  spectra  were  studied  in  order 
to  derive  three  library  spectra.  All  of  the  high  resolution 
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laboratory  or  in  the  field.  In  order  to  make  this  a 
manageable  problem,  some  limitations  and  assumptions  had  to 
be  incorporated. 


First  of  all,  the  source  spectra  are  transported  through 
air  only,  and  not  through  such  materials  as  might  be  found  in 
a  real  transport  problem.  Secondly,  the  detector  counting 
times  are  determined  so  as  to  give  a  predetermined  spectral 
intensity,  whereas  in  reality,  this  flexibility  may  not  exist. 
Thirdly,  the  theoretical  detector  has  a  collecting  area  of  one 
square  centimeter,  and  is  equally  sensitive  to  all  energy 
groups.  Finally,  it  must  always  be  remembered  that  the  Finite 
Elements  Method  itself  is  only  an  approximation  of  the 
transport  solution,  not  an  exact  (analytic)  solution. 

The  assumptions  made  in  this  problem  are  outlined  in 
detail  in  chapter  3,  which  discusses  the  transport  phenomena. 

LAYOUT  OF  CHAPTERS 

The  chapters  are  organized  in  the  order  of  the  subjects 
presented  in  this  introduction.  Chapter  2  begins  with  the 
selection  of  the  sources,  chapter  3  covers  the  transport  of 
the  sources,  and  chapter  4  completes  the  solution  with  a 
detailed  presentation  of  Bayes'  theorem  and  the  counting 
statistics.  Chapter  5  discusses  the  results  and  compares  them 
with  the  method  of  photopeak  identification,  and  the  last 
chapter  outlines  the  conclusions  and  recommendations, 
including  a  discussion  of  the  applications  of  this  method  to 
'real  world'  scenarios. 
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of  the  ideal  detector.  Initially,  the  Poisson  distribution 
will  be  used,  under  the  assumption  that  all  channels  of  the 
detector  collect  counts  independently  of  each  other 
(statistically  independent).  In  order  to  test  the 
sensitivity  of  the  analysis,  a  multinomial  distribution  will 
be  used  instead  of  the  Poisson  distribution,  this  time  under 
the  assumption  that  the  detector  channels  are  statistically 
interdependent. 

Finally,  the  fourth  part  of  this  thesis  problem  entails 
writing  a  subroutine  which  uses  Bayes'  theorem  to  examine  the 
transported  spectra  at  any  given  distance  from  the  source  and 
output  the  probability  that  a  particular  spectrum  is  that  of 
"Source  A"  given  the  transported  spectrum.  This  is  repeated 
for  all  sources  (Source  B,  etc.)  for  the  measured  spectrum  at 
a  number  of  ranges  and  the  results  plotted.  The  Bayesian 
source  liklihoods  can  be  determined  at  ranges  beyond  the 
acceptable  limits  (to  be  defined  in  chapter  V)  of  photopeak 
identification,  and  the  utility  of  the  Bayes'  theorem  method 
will  be  assessed.  The  computer  code(s)  for  the  third  and 
fourth  parts  of  this  thesis  will  be  combined  as  one  and 
written  in  the  BASIC  computer  language  and  run  on  an  Apple  II 
series  computer. 

Limitations 

This  problem  is  designed  to  verify  a  theoretical 
technique  using  numerical  methods.  As  such,  the  thesis  does 
not  pretend  to  duplicate  all  the 


constraints  found  in  a 
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photopeak  analysis  fails.  In  so  doing,  it  will  be  shown  that 
information  about  the  source  can  be  gained  from  the 
downscattered  photons. 
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Approach 

The  problem  will  be  broken  down  into  four  main  parts. 
First,  a  set  of  realistic  ’generic*  source  spectra  will  be 
developed  based  on  the  gamma  ray  spectra  of  isotopes  likely 
to  be  found  in  a  nuclear  source.  A  number  of  different 
sources  will  be  developed  in  order  to  test  the  ability  of 
Bayes*  probability  theory  to  discriminate  between  different 
nuclear  sources. 

The  second  part  of  the  problem  involves  developing  a 
computer  program  to  transport  each  individual  source  spectrum 
for  any  distance  through  the  atmosphere.  The  first  thesis 
done  on  this  topic  attempted  to  transport  gamma  ray  spectra 
using  Monte  Carlo  techniques  (7).  This  method  proved  to  be 
far  too  cumbersome  for  the  given  problem  and  will  not  be 
investigated  further.  Instead,  this  thesis  will  use 
multigroup  diffusion  theory,  which  should  be  a  more 
applicable  transport  technique  (8).  The  multi-group 
diffusion  code  developed  here  uses  the  Finite  Elements  Method 
with  a  one-dimensional  spherical  coordinate  geometry  and  will 
be  written  in  Fortran  77  on  the  AFIT/AD  Harris  800 
minicomputer . 

The  third  major  effort  will  be  to  incorporate  a 
subroutine  to  analyze  the  counting  statistics  at  the  location 
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where  F  is  now  the  function  that  minimizes  the  penalty 
function  P  (3).  Next,  F  is  replaced  by  F  +  6F,  where  6F  is 
arbitrarily  small,  and  equation  (3)  becomes 


P ( F  +  oF)  =  \ 


dV{ DV( F  +  <5F)  »V(F  +  SF) 


(4) 


+  £a(F  +  6F) 2  -  2 ( F  +  6F)S } 


By  multiplying  through  and  collecting  terms  on  the  right  hand 
side,  equation  (4)  becomes 


P(F  +  6F)  =  P(F)  + 


dV{  DV6F  »7F  +  2?F<5F 

-  6FS}  +  0(52F) 


(5) 


The  terms  of  order  52F  in  equation  (5)  can  be  neglected,  as 
they  are  vanishingly  small,  so  that  the  equation  is  of  the 
form 


P(F  +  6F)  =  P(F)  +  P(6F) 

6P  =  P(F  +  <$F)  -  P(F) 


(6) 

(7) 


In  order  to  obtain  the  smallest  penalty,  the  right  hand  side 
of  equation  (7)  must  go  to  zero.  From  equations  (5)  and  (7) 
then 


6P 


dV{Dv6F*VF  +  £aF 5F  -  6FS } 


(8) 


By  using  integration  by  parts  and  Stoke's  theorem,  the 


incremental  penalty  function  becomes 


5P 


fdV6F{-V*(DVF)  +  LaF  -  S} 


(9) 


This  results  in  a  final  form  of  6P,  given  by 


6P 


da_5FDV_F  + 


dV6F{-V*(DVF)  +  IaF  -  S} 


(10) 


Therefore,  to  minimize  this  quantity  (6P),  both  of  the 
integrands  (one  surface  and  one  volume)  must  go  to  zero. 


VF  =  0 

-V«(DVF)  +  ZaF  -  S  =  0 


(ID 

(12) 


It  is  evident  from  the  results  presented  in  equations  (11) 
and  (12)  that  the  penalty  function  given  by  equation  (3) 
accurately  describes  the  diffusion  problem  (comparing 
equations  (1)  and  (12)).  Now  the  problem  has  been  reduced  to 
one  of  minimizing  the  penalty  function  in  equation  (3)  by 
solving  for  the  function  F. 

Digitizing  the  Problem. 

The  methodology  discussed  above  must  now  be  worked  into 
a  form  conducive  to  computer  implementation.  The  finite 
element  nature  of  the  problem  enters  when  the  function  F  is 
represented  by  a  homogeneous  cubic  polynomial  function.  In 
one  dimension,  this  function  is  given  by  equation  (13)  (top 
of  next  page): 
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F(1  1  )  =.Iz  f  (1  1  ) 

l  l  i=iii  i  l 


The  values  of  z  are  the  fluxes  and  currents  (flux 
derivatives)  at  the  mesh  boundaries,  and  hence  represent  the 


ultimate  form  of  the  solution  of  the  finite  element  problem. 
The  functions  f.  in  equation  (13)  are  expressed  in  terms  of 
the  natural  coordinates  1 ^  and  1 ^ • 


1  1  +  1 2  ~  L 


(U) 


In  addition,  1^  has  a  value  of  1  at  the  left  mesh  node  and  0 
at  the  right  node,  whereas  1 ^  is  0  on  the  left  and  one  on  the 
right.  The  mesh  configuration  is  illustrated  in  figure  2 
below.  ^  1  i  ^ 


Figure  2.  One-Dimensional  Version  of  FCl^.lj) 

The  function  F  in  equation  (3)  has  now  been  approximated  in 

terms  of  the  mesh  quantities  F^,  j  ^ ,  F2F  j  2 »  t^ie  fluxes  and 

currents  at  the  nodal  points.  In  its  expanded  form  equation 

(13)  becomes 


F  -  Fi£1d1»12)  +  J1f2<'1l  ,12')  +  F2£3<1r12) 


_\.  V  v*  V 


The  cubic  polynomial  functions  are  given  by  the  following 
equations 
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(19) 


Equation  (3)  can  now  be  expressed  in  terms  of  the  function  F 
in  equation  (15).  Both  the  first  and  second  terms  in 
equation  (3)  can  be  replaced  with  the  aid  of  equations  (16) 
(19).  For  instance,  the  second  term  in  equation  (3)  becomes 


^JdlZa{F1f1  +  j1f2  +  F2f3  +  j2f4>2  (2°) 

Equation  (20)  can  be  expanded  and  the  integral  evaluated  by 
using  the  integral  formula  given  in  the  equation  below  (15). 


-  (p+q  +  1  )  I  <21> 

The  first  term  in  equation  (3)  is  done  in  a  similar  fashion, 
with  the  additional  aid  of  the  relationships 

dl x/dx  =  —  1  / h  (22) 

dl  /dx  =  1/h  (23) 


The  third  term  in  equation  (3)  is  approximated  using 
Simpson’s  intergral  approximation,  since  the  value  of  the 
source  function  S  can  be  calculated  throughout  the  region  of 
interest  ( 5) . 

Once  the  three  terms  in  equation  (3)  have  been 
approximated,  it  is  a  simple  matter  to  get  it  into  a  form 
which  allows  for  computer  coding.  All  three  integral  terms 
have  been  put  in  a  form  allowing  for  evaluation  over  the 
interval  of  interest,  and  equation  (3)  becomes 

P  =  iZ.M.z  -  Z.S  -  X(c .z  -  d)  (24) 

The  first  term  is  simply  a  combination  of  the  first  two 
integral  terms  in  equation  (3),  having  been  expanded  and 
combined  into  a  single  matrix  M.  This  global  matrix  is  a 
combination  of  all  the  local  matrices,  there  being  one  local 
matrix  for  each  mesh  space  used.  Recall  that  the  local 
matrices  were  obtained  by  evaluating  terms  one  and  two  of 
equation  (3)  using  equations  (16)  -  (20). 

The  S  matrix  consists  of  the  values  of  the  third 
integral  term  in  equation  (3),  which  was  evaluated  using 
Simpson's  Rule.  The  last  term  in  equation  (24)  was  added  to 
account  for  the  boundary  conditions  of  the  problem. 

Finally,  the  z  matrix  is  an  array  containing  the  values 
of  the  flux  and  current  at  each  node.  Hence,  the  problem 
remains  to  solve  for  the  values  of  the  z  matrix  by  minimizing 


equation  (24).  By  taking  the  partial  derivative  of  P  in 
equation  (24)  with  respect  to  z  and  setting  it  equal  to  zero 
as  well  as  the  partial  derivative  of  M  with  respect  to  X ,  a 
series  of  linear  equations  is  obtained.  The  equations, 
expressed  in  matrix  form  are  rearranged  and  ultimately  used 
to  solve  for  the  z  matrix 

z  =  M~lS  +  M"1c1X1  +  M_1c2X2  (25) 

The  transport  code  is  written  to  solve  just  this  set  of 
equations . 

The  One-Dimensional,  Single  Group  Problem. 

As  a  first  step  to  running  the  transport  code  using  the 
Finite  Element  Method,  a  simple  case  was  programmed.  This 
case  involves  a  one  dimensional  "slab  geometry"  problem  with 
a  known  source  function  S,  and  a  known  fluence  solution  (12) 
The  diffusion  equation  for  this  geometry  becomes: 

-D(if^Ll)2  +  £aF(  x  )  =  S  ( x )  (26) 


and  the  final  equation  for  the  penalty  function  becomes 

a 

P  =  +  £aF2  "  2FS )  (27) 

b 


This  is  precisely  the  problem  which  was  derived  in  the  last 
section,  using  the  one-dimensional  cubic  polynomials 


(equation  15).  Since  the  fluence  as  a  function  of  distance 
is  already  known,  it  becomes  an  easy  matter  to  check  the 
results  of  the  program. 

The  program  was  tested  on  the  case  where  the  fluence  and 
source  functions  are  given  by 

F ( x )  =  x4  (28) 

S ( x )  =  -1 2x  2  +  x  4  (29) 


The  finite  element  solution  (using  the  source  function)  is 
graphed  along  with  the  actual  solution  in  figure  3. 


Figure  3.  Finite  Elements  Solution  Verification 


As  can  be  seen  in  figure  3,  the  match  provided  by  the  Finite 
Elements  Method  is  very  good  indeed. 

The  Three-Dimensional.  Multi-Group  Problem 


Following  the  code  verification  in  one  dimension,  it  is 
necessary  to  expand  the  code  into  the  three-dimensional 
(spherical  geometry)  multigroup  problem.  Since  assumption 
number  (5)  assumes  radial  symmetry,  the  only  dependence  is 
on  the  variable  r.  As  is  explained  in  appendix  C,  the 
three-dimensional  problem  can  be  placed  into  a  form  which 
allows  the  use  of  the  one-dimensional  approximating 
polynomials,  given  by  equations  (15)  -  (19).  Hence  the 
problem  needs  only  to  be  modified  to  accept  the  multiple 
energy  groups. 

The  transition  from  single  to  multigroup  is  very 
straight-forward.  The  program  is  modified  to  calculate  the 
solutions  for  flux  and  current  on  the  nodes  (the  z  matrix) 
one  time  for  each  group.  This  involves  setting  up  the  global 
matrix  once  for  each  group  as  well  as  modifying  the  source 
term  used  in  the  Simpson's  approximation.  Beginning  with  the 
highest  energy  group  (group  1),  the  diffuse  fluence  solution 
is  due  only  to  inscattering  of  the  uncollided  photons.  As 
the  program  solves  each  consecutively  lower  energy  group,  the 
source  term  is  modified  to  account  for  the  downscat tered 
fluence  from  higher  energy  groups  as  well  as  the  inscattered 
contribution.  The  iterative  process  proceeds  in  the 
direction  of  decreasing  energy  groups,  and  hence  upscattering 
from  lower  energy  groups  to  higher  energy  groups  is  not 
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allowed . 


In  order  to  verify  the  full  scale  radiation  transport 
code,  a  known  function  similar  to  equations  (28)  and  (29)  was 
input,  only  in  this  case  it  was  a  simple  exponential 
function.  Again,  the  code  output  values  are  very  close  to 
the  actual  solution.  A  more  important  bench  mark  is  attained 
by  solving  the  diffusion  equation  analytically  for  a  simple 
case.  The  analytic  solution  is  outlined  in  appendix  D,  along 
with  the  value  returned  by  the  transport  code.  As  can  be  seen 
in  appendix  D,  the  answer  is  acceptable,  given  only  a  one 
mesh  space  approximation. 

Source  Transport 

The  next  step  in  the  process  is  the  actual  transport  of 
each  library  source  to  the  desired  ranges.  Each  of  the 
sources  was  reconstructed  at  the  desired  ranges,  yielding  a 
degraded  set  of  source  spectra  due  to  scattering  and 
absorption  losses.  The  two  types  of  transported  sources  are 
discussed  below. 

The  Measured  Spectra 

In  order  to  simulate  measured  spectra  at  various  ranges 
from  the  sources,  it  is  necessary  to  transport  the  library 
spectra  in  two  different  fashions.  The  first  method  involves 
using  the  transport  code  (developed  above)  alone  to  determine 
the  relative  intensities  of  each  library  spectrum  at  a  number 
of  different  ranges.  These  spectra  represent  the  ideal 
measured  spectra,  since  the  transport  code  does  not  introduce 


25 


any  random  processes  or  counting  statistics  by  itself.  These 
"reference  spectra"  thus  represent  good  data  sets  which  would 
be  obtained  from  an  enormous  number  of  measurements  in  the 
field  (i.e.,  the  true  mean  in  each  channel). 

The  second  method  involves  taking  the  reference  (ideal) 
spectra  mentioned  above  and  converting  each  one  into  a 
measured  spectrum  by  introducing  random  counting  statistics. 
The  statistics  subprogram  introduces  random  variations  in  the 
counting  data  in  li^e  with  the  appropriate  statistical 
distribution.  Two  different  distributions  will  be  employed  in 
deriving  the  measured  source  spectra:  the  Poisson 
distribution  and  the  multinomial  distribution.  Since  the 
measured  sources  are  actually  generated  during  the  analysis 
part  of  the  problem,  the  methodology  will  be  covered  in 
detail  in  chapter  IV  and  its  appendices. 

There  are  now  two  different  sets  of  spectra  which  will 
be  used  in  the  analysis  process.  The  first  set  consists  of  a 
good  group  of  measured  spectra;  the  "reference  spectra," 
which  was  obtained  by  the  transport  of  the  library  spectra. 
These  reference  spectra  are  equivalent  to  a  set  of  spectra 
determined  by  many  different  measurements  of  the  sources. 

The  reference  sources  have  been  plotted  out  at  ranges  of  50, 
100,  200,  300,  400  and  500  meters,  and  are  presented  in 
figures  4  through  9.  It  is  important  to  note  that  the 
reference  source  spectra  are  in  47rr2  fluence  units,  so  they 
can  be  easily  compared  with  the  library  sources.  The  actual 
values  of  fluence  for  the  reference  sources  are  tabulated  in 


tables  III-l  through  III-6. 

The  second  set  of  source  spectra  represents  a  less  than 
ideal  group  of  spectra;  the  "measured  spectra,"  obtained  by 
the  transport  of  the  library  spectra  with  the  addition  of 
counting  statistics  considerations.  This  group  is  equivalent 
to  a  set  of  spectra  obtained  by  taking  a  single  measurement 
of  the  sources.  Since  the  measured  sources  are  generated 
during  the  analysis  process,  they  will  be  listed  in  chapter 
IV. 

Background  Considerations 

In  addition  to  the  reference  spectra  and  measured 
spectra  discussed  above,  the  library  spectra  were  generated  a 
third  time,  but  this  time  background  counts  were  accounted 
for.  The  background  spectrum  provided  with  the  high 
resolution  data  was  used  to  determine  the  mean  number  of 
background  counts  present  in  each  channel  of  the  reference 
source  spectra.  The  measured  source  spectra  are  generated  in 
a  manner  similar  to  the  "no  background  case,"  but  now  a 
separate  subroutine  is  used  to  generate  random  background 
counts  in  line  with  the  appropriate  statistics.  These 
background  counts  are  then  added,  channel  by  channel,  to  the 
measured  source  spectrum  of  interest.  Again,  the  generation 
of  the  measured  source  spectra  is  discussed  at  length  in 
chapter  IV. 
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Figure  4.  Reference  Source  A  at  50,  100  &  200  meters 
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Figure  6.  Reference  Source  B  at  50,  100  &  200  meters 


Table  I V— 1 :  Sample  Prior  Distribution 


Source 


Number  in  Sample 


PCAi) 


1  5  0.5 

2  2  0.2 

3  3  0.3 


The  second  term  in  the  numerator  of  equation  (30), 
P(X|Ai),  is  the  probability  of  obtaining  the  measured 
spectrum  X  given  the  reference  spectrum  Ai.  This  term  is  a 
calculated  quantity,  and  will  vary  depending  on  the  type  of 
counting  statistics  used. 

Poisson  Distribution. 

The  first  method  of  determing  the  values  of  P(X|Ai)  will 
assume  that  the  counts  are  collected  in  a  Poisson 
distribution.  The  assumption  then  is  that  each  one  of  the 
eighteen  channels  in  the  low  resolution  detector  is 
statistically  independent  from  the  others.  The  Poisson 
distribution  is  a  particularly  valid  one  in  counting  problems 
when  the  average  number  of  counts  is  much  smaller  than  the 
total  possible  number  (i.e.,  smaller  than  the  number  of  gamma 
photons  emitted  by  the  entire  nuclear  source)  (4).  Counting 
data  collected  in  a  Poisson  distribution  is  obtained  by 
running  the  detector  for  an  arbitrary  length  of  time,  which 
generally  yields  a  different  number  of  total  counts  each 
time.  This  is  in  contrast  to  the  multinomial  distribution, 
which  will  be  covered  later  in  this  chapter. 


represents  the  prior  distribution  of  the  three  different 
sources.  As  a  starting  point,  the  three  values  will  be 
assumed  equal,  with 


P(A1)  =  0.3333 

(33) 

P(A2)  =  0.3333 

(34) 

P(A3)  =  0.3333 

(35) 

In  other  words,  it  is  assumed  that  each  of  the  three  sources 
occurs  equally  often.  This  quantity  can  be  varied  to  test 
the  sensitivity  of  the  posterior  distribution.  In  the  case 
of  actual  source  measurements  in  the  field,  the  prior 
distribution  would  presumably  be  know. 

As  a  simple  example  of  the  prior  distribution,  assume 
ten  containers  of  nuclear  materials  have  had  the  identifying 
labels  removed.  Further,  assume  that  the  only  way  to 
determine  the  identity  of  the  nuclear  material  in  each  of  the 
containers  is  to  measure  the  spectrum  of  each  with  a 
multichannel  gamma-ray  detector  at  a  range  of  100  meters.  In 
order  to  use  Bayes'  Theorem,  the  prior  distribution  must  be 
determined.  In  this  example  problem,  the  total  number  of 
each  type  of  material  is  known,  so  it  is  a  simple  matter  to 
determine  the  prior  distribution,  as  seen  in  table  IV-1.  In 
the  main  thesis  problem,  it  is  assumed  that  the  prior 
distribution  is  also  a  given. 
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which  is  the  joint  probability  of  the  variables  X  and  Ai 
(17).  Or,  in  other  words 


P(A  ,X)  =  P  (  X  ,  A  )  (32) 

The  denominator  in  equation  (30)  is  the  sum  of  all  the 
combined  likelihood  functions  and  prior  distributions.  The 
summation  is  over  the  total  number  of  reference  sources, 
which  is  three  in  this  problem. 

As  can  be  seen  in  equation  (30),  Bayes'  theorem  allows 
for  the  combination  of  prior  information  (the  prior 
distribution  of  sources)  and  the  measured  data  (the 
likelihood  functions)  to  determine  the  posterior 
distribution,  given  by  P(Ai|X).  The  posterior  distribution 
gives  the  probability  (likelihood)  of  having  source  Ai, 
given  the  measured  spectrum  X.  For  a  given  measured  source 
X,  P ( A 1 | X )  is  the  probability  that  the  source  AI  is  actually 
the  one  being  measured.  The  power  of  this  method  is  that  it 
clearly  allows  for  the  combination  of  all  available 
information  in  arriving  at  the  posterior  distribution. 

Applied  Form 

Having  reviewed  Bayes'  Theorem  in  its  general  form,  the 
next  task  is  to  examine  it  specifically  in  relation  to  this 
problem.  The  equation  will  be  examined  term  by  term,  and  the 
specific  quantities  will  be  identified  as  they  apply  to  this 
problem . 

As  was  mentioned  above,  the  P(Ai)  term  in  equation  (30) 


(30) 


P(X|A.)  P(A.) 

I  P(X|A  )P(A 
j=1  J  J  J 

The  variable  Ai  (where  i  =  1,3)  represents  the  known  data,  or 
in  this  case,  the  reference  sources,  and  X  represents  the 
given  sample  data,  or  the  measured  source.  Therefore,  P(Ai) 
is  the  probability  of  occurence  of  the  reference  source  Ai. 
Quite  simply,  P(Ai)  represents  the  proportion  of  occurrence 
of  source  Ai,  so  that  if  there  exists  a  group  of  ten  nuclear 
sources,  of  which  four  are  known  to  be  of  type  AI  and  six  of 
type  A2,  then  P(A1)  =  0.40  and  P(A2)  =  0.60.  This 
probability  distribution  is  defined  as  the  "prior 
distribution."  The  prior  distribution  is  the  unconditional 
probability  of  Ai  occurring. 

Given  that  X  represents  the  observed  or  measured  value 
of  a  source,  then  the  P ( X | A i )  term  in  equation  (30)  is  the 
probability  that  the  measured  source  spectrum  X  is  due  to  the 
reference  spectrum  Ai.  Given  the  source  spectrum  AI  then, 
P(X| AI)  is  the  probability  that  a  measured  spectrum  X  is 
caused  by  the  spectrum  AI.  This  probability  is  defined  as 
the  likelihood  function,  and  is  a  discrete  value  in  this 
problem,  since  there  is  a  predefined  number  of  sources.  Note 
that  the  numerator  in  equation  (30)  can  also  be  expressed  as 


P(A.|X) 


P(A  ,X)  =  P ( X | A  )  P( A  .  ) 


(31) 


IV. 


Identification  and  Analysis  of  Source  Spectra 


In  chapter  II,  the  library  source  spectra  were  derived, 
and  the  methodology  for  creating  the  reference  spectral 
sources  and  the  measured  spectral  sources  was  described.  In 
chapter  III,  the  transport  code  was  developed  using  the 
Finite  Element  Method  (FEM).  This  code  was  employed  to 
generate  the  reference  source  spectra  discussed  in  chapter 
III.  Now  that  all  the  raw  data  is  available,  a  technique  must 
be  defined  by  which  the  measured  sources  can  be  identified 
and  analyzed.  The  methods  of  identification  and  analysis 
will  be  discussed  in  this  chapter. 

Bayes  Theorem 

General  Form. 

The  method  of  identifying  sources  in  this  thesis  will  use 
Bay es'  Theorem .  In  the  general  sense,  the  term  Bayesian  is 
used  to  describe  an  approach  to  combining  information  that 
uses  the  sample  information  (i.e.,  data  collected)  as  well  as 
other  available  information  (17).  In  the  case  of  spectral 
source  identification,  the  sample  information  consists  of  the 
measured  spectra,  and  the  other  information  is  the  set  of 
reference  spectra. 

Bayes'  theorem  allows  one  to  make  an  assessment  of  the 
likelihoods  of  the  occurence  of  certain  events,  given  the 
sample  information.  In  terms  of  the  variables  X  and  A  (top 
of  next  page): 


TABLE  III-6 

Relative  Intensity  Distribution  of  Reference  Source  C 

(counts  /  photon) 


up  Number 


300  meters 


400  meters 


500  meters 


1 

0.0 

0.0 

0.0 

2 

0.0 

0.0 

0.0 

3 

0.0 

0.0 

0.0 

4 

0.0 

0.0 

0.0 

5 

0.0 

0.0 

0.0 

6 

2.78E-15 

9.35E-16 

3.58E-16 

7 

1.45E-14 

4.73E-15 

1.81E-15 

8 

2 . 28E-14 

7 . 23E-15 

2.84E-15 

9 

1.02E-14 

3.96E-15 

2.29E-15 

0 

3.37E-14 

1.07E-14 

4.66E-15 

1 

1.75E-14 

6.23E-15 

3.43E-15 

2 

3.54E-14 

1.17E-14 

5.81E-15 

3 

1 . 29E-13 

4.04E-14 

1 .90E-14 

4 

1.52E-13 

4. 5 IE- 14 

2. 18E-14 

5 

1.14E-12 

3.36E-.3 

1.54E-13 

6 

1.69E-11 

6.31E-12 

3.15E-12 

7 

4.00E-11 

1.62E-11 

9.24E-12 

.8 

6.75E-12 

2.73E-12 

1 . 58E-12 

^  . 
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TABLE  III-5 


Relative  Intensity  Distribution  of  Reference  Source  C 

(counts  /  photon) 


GrouD  Number 


50  meters 


100  meters 


200  meters 


6 

3.61E-13 

6.98E-14 

1 .04E-14 

7 

2.02E-12 

3.85E-13 

5 . 61E-14 

8 

3.40E-12 

6.43E-13 

9.14E-14 

9 

2.37E-13 

1.08E-13 

3.05E-14 

0 

5.23E-12 

9.76E-13 

1.37E-13 

1 

5.42E-13 

2.31E-13 

5.81E-14 

2 

4.76E-12 

9.22E-13 

1.39E-13 

3 

2.67E-11 

4.55E-12 

5.60E-13 

4 

2.76E-11 

5.08E-12 

6.75E-13 

5 

3.56E-10 

5.55E-11 

5.71E-12 

6 

1.29E-9 

3.00E-10 

5.67E-11 

7 

1.31E-9 

4.49E-10 

1 .  17E-10 

8 

3.41E-10 

8.99E-11 

2.02E-11 

TABLE  III-4 

Relative  Intensity  Distribution  of  Reference  Source  B 

(counts  /  photon) 


Group  Number 


300  meters 


400  meters 


500  meters 


4.49E-15 

1 .51E-15 

5.80E-16 

1.60E-14 

5.29E-15 

2.06E-15 

1.15E-14 

3.91E-15 

1.77E-15 

1.04E-I4 

3.76E-15 

I.95E-15 

1.30E-14 

4.85E-15 

2.65E-15 

I.33E-14 

4.45E-15 

2.23E-I5 

9.46E-14 

2.39E-14 

8.04E-15 

5.84E-12 

1.53E-13 

5.60E-14 

1.14E-12 

2.71E-I3 

1.01E-13 

3.97E-12 

1.59E-12 

9.56E-13 

2.45E-11 

1.05E-11 

6.19E-12 

4.54E-11 

2.08E-11 

1 .30E-11 

7.53E-12 

3.46E-12 
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TABLE  I I 1-2 

Relative  Intensity  Distribution  of  Reference  Source  A 

(counts  /  photon) 


Number 


300  meters 


400  meters 


500  meters 


TABLE  III-l 

Relative  Intensity  Distribution  of  Reference  Source  A 

(counts  /  photon) 


Group  Number 


50  meters 


100  meters 


200  meters 


0.0  0.0  0.0 

0.0  0.0  0.0 

0.0  0.0  0.0 

0.0  0.0  0.0 

0.0  0.0  0.0 

5.41E-13 

1.05E-13 

1.57E-14 

1.92E-12 

3.69E-13 

5.42E-14 

4.02E-12 

7.56E-13 

1.07E-13 

5.04E-12 

9.70E-13 

1.42E-13 

1.60E-11 

2.92E-12 

3.91E-13 

2.23E-11 

4.03E-12 

5.15E-13 

5.37E-11 

9.51E-12 

1 . 19E-12 

1.05E-9 

1 .68E-10 

1 .76E-11 

7.03E-10 

1 .37E-10 

1 .85E-11 

2.40E-10 

8.97E-11 

2 . 67E-1 1 

5. 14E-I0 

2.35E-10 

8 . 41E-1 1 

3.56E-10 

2. 15E-10 

1 .03E-10 

1 .93E-10 

4.87E-11 

1.72E-11 
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R  =  500  METERS 
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Figure  9.  Reference  Source  C  at  300,  400  & 
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REFERENCE  SOURCE  B 
R  =  300  METERS 


REFERENCE  SOURCE  B 
R  =  500  METERS 


The  Poisson  distribution  is  given  by 


& 


x 

P( x , a )  =  exp(-a) 


(36) 


where  the  variable  a  is  the  mean  number  of  counts  occurring 
in  the  channel  of  interest,  and  x  is  the  observed  number  of 
counts  in  the  channel  (4).  Or,  in  the  case  of  a  multichannel 
detector,  equation  (36)  becomes 


p(xk'“ik)  - 


<aik)*k 


exp(-aik) 


(37) 


where  k  denotes  the  channel  number,  which  ranges  from  one  to 
eighteen.  In  order  to  get  an  expression  for  P(X|Ai),  the 
individual  channel  probabilities  must  be  combined  into  a 
single  probability.  This  is  achieved  by  noting  that  since 
the  channels  are  statistically  independent  of  one  another, 
the  probabilities  multiply.  Hence  the  relationship  becomes 


18  (a  . ,  )Xk 

P(X|Ai)  =  II  XX ,  exp(-aik)  (38) 

k=1  k  ' 


The  mean  number  of  counts  in  channel  k  is  given  by  aik,  and 
comes  from  the  reference  spectrum,  discussed  in  chapter  III 
(ref.  figures  4  -  9).  The  final  form  of  equation  (38)  then 
becomes  (top  of  next  page): 
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18  (S  +  b  )xk 

P(X|A  )  =  n  —  -  exp(-(Sik  +  Bk))  (39) 

k=1  k  ’ 

Here  Sik  represents  the  mean  number  of  counts  from  reference 
spectrum  i  (i  =  1,3)  in  channel  number  k  (k  =  1,18)  and  Bk  is 
the  mean  number  of  background  counts  in  channel  k.  Recall 
that  initially,  the  analysis  is  done  with  zero  background 
counts  in  all  channels  of  all  sources.  Xk  in  equation  (39) 
represents  the  number  of  counts  in  channel  k  observed  in  the 
measured  spectrum. 

The  Poisson  version  of  the  Bayes'  theorem  analysis  was 
coded  in  the  BASIC  computer  language  and  run  on  an  Apple  II 
series  computer.  The  program  listing  and  description  is 
included  in  appendix  H. 

The  value  of  Xk,  obtained  from  the  measured  source 
spectrum,  is  obtained  using  a  random  number  generator  and  the 
Poisson  distribution.  Details  of  this  process  are  covered  in 
appendix  G.  On  the  other  hand,  the  values  of  Sik  come  from 
the  reference  source  spectra,  which  were  discussed  in  chapter 
III.  The  method  of  obtaining  numerical  values  for  the 
reference  sources  is  discussed  in  appendix  F. 

Gaussian  Distribution. 

In  the  limit  of  a  large  number  of  counts,  the  Poisson 
distribution  approaches  the  normal  or  Gaussian  distribution, 
which  is  given  by  (top  of  next  page): 
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(40) 


P(X|A.) 


18 

n 

k=1 


exp(-Hxk  -  (S.k  +  Bk))/0.k) 

(2Tra2ik)  ^ 


Here,  oik  is  the  standard  deviation  of  the  counts  of  source  i 
in  channel  k.  The  standard  deviation  squared  is  therefore 
equal  to  the  number  of  counts  observed  in  the  channel  of 
interest,  Xk.  All  other  terms  in  equation  (40)  are  defined 
the  same  as  in  the  Poisson  distribution. 

In  practice  (in  this  problem)  the  Gaussian  distribution 
is  used  anytime  the  counts  measured  in  any  one  channel  exceed 
34,  since  the  factorial  in  the  denominator  of  the  Poisson 
distribution  (equation  (39))  exceeds  the  overflow  limits  of 
the  computer  (1E+38).  The  Gaussian  version  of  the  Bayes' 
theorem  analysis  was  also  coded  in  BASIC,  and  is  included  in 
appendix  H. 

A  new  measured  source  spectrum  is  generated  each  time 
the  analysis  code  is  run,  and  hence  determines  the  values  for 
Xk  in  equation  (40).  The  method  of  obtainining  a  measured 
source  spectrum  with  a  Gaussian  distribution  is  contained  in 
appendix  G.  The  method  of  obtaining  reference  source  spectra 
is  the  same  as  in  the  Poisson  version. 

Multinomial  Distribution. 

The  second  main  method  of  performing  the  Bayes'  theorem 
analysis  employs  the  multinomial  distribution.  Unlike  the 
Poisson  (and  Gaussian)  distribution,  the  counts  in  the 


different  channels  of  a  multinomial  distribution  are  all 


interdependent  (2).  To  obtain  a  multinomial  distribution  in 
a  multichannel  analyzer,  counts  are  collected  until  a  certain 
predetermined  total  is  reached,  at  which  time  the  detector  is 
turned  off.  By  collecting  in  this  manner,  the  number  of 
counts  collected  in  one  channel  becomes  dependent  on  the 
counts  in  all  other  channels  (10).  This  dependence  follows 
quite  simply  because  each  count  collected  in  one  channel 
decreases  by  one  the  total  counts  available  for  all  other 
channels.  In  terms  of  equation  (30)  variables,  the 
multinomial  distribution  is  given  as 

P(X|A  )  =  n  ( f  ..  )Xk  (41) 

Ji  xk!  k=1 

Here  Xk  again  represents  the  measured  number  of  counts  in 
channel  k,  subject  to  the  constraint 

1  8 

I  \  -  N  (42) 

k=1 

where  N  is  the  the  total  number  of  counts  in  all  channels 
combined  (13).  The  variable  fik  is  the  fraction  of  the 
counts  of  source  i  occurring  in  channel  k.  The  fik's  are 
obtained  from  the  reference  source  spectra  and  the 
methodology  is  discussed  at  length  in  appendix  F.  The  fik 
values  are  subject  to  the  constraint 


for  each  source. 

As  was  true  in  the  Poisson  and  Gaussian  analysis  codes, 
the  multinomial  analysis  code  is  contained  in  appendix  H. 

The  method  of  obtaining  the  measured  source  (Xk)  is  discussed 
in  appendix  G,  and  the  reference  sources  (fik)  in  appendix  F. 


V .  Result s  a nd  Discussion 


Having  the  reference  sources  in  hand,  the  problem 
remains  to  analyze  the  measured  sources  at  the  ranges  of  50, 
100,  200,  300,  400  and  500  meters  from  the  source.  The 
sections  in  this  chapter  will  be  presented  in  order  of  the 
analysis  runs.  First  the  Bayes*  theorem  analysis  is  done 
using  Poisson  statistics  with  no  background  contribution. 
Secondly,  the  multinomial  version  of  Bayes'  theorem  is  used, 
again  with  no  background.  Thirdly,  the  Poisson  version  will 
be  run  again,  but  this  time  using  the  background  contribution 
discussed  earlier.  Finally,  an  unknown  source  will  be 
"measured"  in  one  pass  at  a  distance  of  50  meters,  to  see  how 
the  Bayes'  posterior  distribution  behaves  with  a  nonlibrary 
source  spectrum. 


Bayes'  Theorem  With  Poisson:  No  Background 

The  Bayes  analysis  code  using  the  Poisson  statistics 
package  was  run  at  all  ranges  using  each  of  the  three  sources 
as  a  measured  source  once  at  each  range.  Hence,  this 
constituted  18  runs  using  the  Poisson  and  Gaussian  codes 
(depending  on  the  level  of  counts,  ref.  chapter  IV).  The 
baseline  level  of  counts  was  assigned  as  1000  counts  at  50 
meters.  This  corresponds  to  a  collecting  time  of 
approximately  30  minutes,  using  the  one  square  centimeter 
detector  assumed  in  this  problem.  An  equal  counting  time  was 
assumed  for  each  measurement,  in  line  with  Poisson  counting 


50 


statistics.  The  total  number  of  reference  counts  in  each 


source  and  the  measured  source  values  for  each  run  were 
determined  as  discussed  in  appendices  F  and  G. 

The  results  of  the  posterior  distribution  are  presented 
in  table  V-l.  As  can  be  seen  in  the  posterior  distribution 
for  measured  source  A,  the  probability  of  positively 
identifying  the  measured  source  drops  below  the  99.7%  level 
near  the  200  meter  mark,  but  stays  above  70%  all  the  way  out 
to  500  meters.  Measured  sources  B  and  C  do  not  stay  above 
the  70%  level,  although  source  C  remains  well  above  the  90% 
level  out  to  400  meters.  Part  of  the  reason  for  the 
consistantly  high  values  with  source  A  is  because  of  the 
relative  strength  of  the  source.  As  can  be  seen  in  table 
V-3,  measured  source  A  has  a  total  of  20  counts  collected, 
which  is  twice  the  number  of  counts  in  both  sources  B  and  C 
at  500  meters.  This  follows  because  library  source  A  has 
more  high  energy  features  relative  to  sources  B  and  C,  and 
high  energy  photons  are  degraded  less  severely  than  are  low 
energy  photons  (reference  cross  section  data  in  appendix  B). 

When  comparing  these  results  to  the  method  of  photopeak 
identification  (reference  appendix  A),  it  is  evident  that  the 
uncollided  photons  are  rapidly  degraded  by  the  exponential 
scattering  term,  and  the  peaks  soon  drop  below  the  3  sigma 
level.  As  can  be  seen  in  appendix  A,  the  characteristic 
photopeak  of  source  B  has  already  degraded  below  the  3  sigma 
level  at  50  meters,  and  the  photopeaks  of  sources  A  and  C 
drop  below  3  sigma  shortly  beyond  50  meters.  Therefore, 
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Table  V-2 


Measured  Source  A 

(Poisson 

Distribution) 

channel 

50  m 

100  m 

200  m 

number 

counts 

counts 

counts 

1 
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0 

0 
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0 

0 

0 
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0 

4 

0 
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0 
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0 
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222 
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5 
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77 
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9 

16 
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73 

26 

17 
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69 

32 

18 

60 
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Table  V -3 


Measured  Source  A  (Poisson  Distribution) 


channel  300  m  400  m 

number  counts  counts 


500  m 
counts 


Table  V-4 


Measured  Source  B  (Poisson  Distribution) 


channel  50  m  100  m 

number  counts  counts 


200  m 
counts 


Table  V-5 

Measured  Source  B  (Poisson  Distribution) 


channel 


300  m 


400  m 


500  m 


number 


counts 


counts 


counts 


Table  V-6 


Measured  Source  C  (Poisson  Distribution) 


channel  50  m  100  m 

numb e  r  counts  counts 


200  m 
counts 


Table  V-7 


Measured  Source  C  (Poisson  Distribution) 

500  m 
counts 


channel  300  m  400  m 

number  counts  counts 


eventually  the  probability  in  equation  (40)  becomes  zero 
(exceeding  the  computer  limits  of  IE-38).  This  leads  to  a 
zero  term  in  the  denominator  of  equation  (30)  of  chapter  IV, 
leaving  an  undefined  Bayes'  posterior  distribution.  This 
problem  is  even  more  likely  to  occur  at  higher  count  levels, 
that  is,  in  regions  close  to  the  source  where  photopeak 
identification  is  at  its  best.  Hence,  when  the  Bayes' 
theorem  method  yields  "undefined"  probabilities,  it  serves  as 
an  indicator  that  an  unknown  source,  that  is,  a  source  not  in 
the  library,  has  been  encountered.  Other  methods  must  then 
be  used  to  determine  the  composition  of  the  new  source. 

Applying  the  Results 

What  can  be  learned  from  these  results  that  can  be 
applied  to  experiments  in  the  field?  First  of  all,  it 
appears  from  the  contrast  between  the  Poisson  and  multinomial 
statistics  that  it  is  far  better  to  run  the  detector  as  long 
as  possible  in  order  to  optimize  the  number  of  counts 
(Poisson).  If  the  detector  is  automatically  shut  off  after 
reaching  a  certain  number  of  counts,  much  information  is  lost 
from  the  stronger  sources.  This  effect  was  observed  in  the 
Poisson  runs  without  background.  Indeed,  even  though  all  of 
the  sources  used  in  this  analysis  were  assumed  to  be  of  equal 
strength,  a  slight  statistical  increase  in  the  number  of 
counts  (source  A)  was  significant  in  terras  of  positive 
identi f ication . 


Another  point  of  interest  uncovered  in  the  analysis  is 


Table  V-12 


Measured  Source  B 

(Poisson  Distribution  with  Background) 


channel 


number 

50  m 

100  m 

200  m 

300  m 

400  m 

500m 

1 

3 

1 

2 

1 

1 

0 

2 

1 

0 

1 

3 

1 

0 

3 

0 

0 

0 

1 

1 

2 

4 

2 

1 

0 

0 

0 

1 

5 

2 

2 

0 

1 

2 

1 

6 

0 

1 

4 

1 

7 

3 

7 

5 

5 

5 

5 

5 

6 

8 

5 

5 

4 

4 

5 

5 

9 

22 

20 

23 

21 

20 

19 

10 

24 

23 

22 

23 

24 

22 

11 

240 

242 

241 

243 

241 

240 

12 

31 

24 

24 

22 

22 

23 

13 

90 

50 

41 

44 

42 

39 

14 

153 

56 

39 

39 

39 

40 

15 

192 

88 

63 

62 

60 

61 

16 

488 

241 

168 

150 

149 

149 

17 

370 

197 

115 

94 

87 

86 

18 

111 

59 

45 

41 

38 

38 
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Table  V-ll 


Measured  Source  A 

(Poisson  Distribution  with  Background) 

channel 

number  50  m  100  m  200  m  300  m  400  m 


PROBABILITY  PROBABILITY  PROBABILITY 


1 

Table  V 

-10 

] 

Bayes*  Posterior 

Distribution 

U 

n 

(Poisson  with 

background  ) 

N 

1.  Measured: 

Source  A 

RANGE 

(m) 

PCAl.pO 

PCA2IX) 

P(A3 | X)  £ 

50 

1.000 

0.000 

0.000 

100 

1.00 

7 . 463E-1 1 

4 . 244E-19 

200 

5.384E-1 

3.026E-1 

1.591E-1  t 

300 

3.783E-1 

3.424E-1 

2.793E-1 

400 

3.642E-1 

3.339E-1 

3.019E-1 

500 

3.611E-1 

3.319E-1 

£ 

3.070E-1  | 

2.  Measured 

Source  B 

RANGE 

Ssl 

P(4_iJ  .XJL 

P( A2 1 X) 

P(  A3  1  X  ) 

50 

0.000 

1.000 

0.000  [ 

100 

8.479E-15 

9.886E-1 

1.140E-2 

200 

1.488E-1 

4.297E-1 

4.215E-1 

300 

2.841E-1 

3.739E-1 

3.420E-1  £ 

400 

3.104E-1 

3.464E-1 

3.432E-1  '> 

'  * 

500 

3.368E-1 

3.380E-1 

3.252E-1  ’.I; 

3.  Measured 

Source  C 

& 

k 

RANGE 

(m) 

P( A 1 1 X ) 

PC  A2 1 X) 

P(  A3  1  X) 

50 

0.000 

0,000 

1.000 

100 

0.00 

1.691E-3 

9.983E-1  f 

4.247E-1 

200 

1 . 571E-1 

4.182E-1 

300 

2.814E-1 

3.447E-1 

3.739E-1 

400 

2.764E-1 

3.514E-1 

3.722E-1  £ 

500 

3.060E-1 

3.442E-1 

3.499E-1  V 

» 
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probabilities  for  positive  source  identification  degrade  much 
more  rapidly  than  in  the  case  with  no  background.  All  three 
sources  maintain  high  probabilities  of  positive 
identification  out  to  100  meters,  at  which  time  the 
background  level  overwhelms  the  number  of  source  counts.  But 
as  can  be  seen  in  tables  V-ll  -  V-13,  the  source  counts  in 
each  channel  are  already  becoming  obscured  by  background  at 
100  meters  (compare  with  tables  V-2  -  V-7).  In  cases  with  a 
large  amount  of  background  then,  Bayes'  theorem  appears  to 
lose  information  quite  rapidly,  but  still  offers  some 
improvement  over  the  method  of  photopeak  identification. 

The  Unidentified  Source 

One  of  the  severe  limitations  of  Bayes'  theorem  applied 
to  spectral  identification  appears  to  be  in  the  realm  of 
identifying  unknown  sources,  or  sources  simply  not  contained 
in  the  source  library.  When  an  unknown  source  was  formulated 
which  differed  considerably  from  the  library  sources,  the 
Poisson  (Gaussian)  analysis  failed  to  return  any  numerical 
information  whatsoever. 

The  problem  is  identified  when  examining  equation  (40) 
in  chapter  IV.  When  using  the  normal  distribution  (i.e.,  in  a 
region  of  reasonable  count  levels),  any  excessive  deviation 
from  the  reference  value  in  a  channel  results  in  an  extremely 
small  value  of  P(X  |Ai),  due  to  the  exponential  terra.  This 
problem  is  compounded  if  every  one  of  the  eighteen  channels 
differs  from  each  channel  of  the  reference  sources,  since 


Table  V-9 


• 

Relative 

Intensity  and  Count 

Levels  of  Background 

Channel 

Counts/ Dhoton 

Counts 

• 

1 

0.0 

0.0 

2 

0.0 

0.0 

• 

3 

0.0 

0.0 

4 

0.0 

0.0 

5 

0.0 

0.0 

6 

2.58E-3 

1.29 

7 

1.01E-2 

5.04 

8 

9.17E-3 

4.59 

# 

9 

4.18E-2 

20.89 

10 

4.56E-2 

22.82 

11 

3.64E-2 

240.00 

• 

12 

4.51E-2 

22.56 

13 

8.11E-2 

40.57 

14 

7.79E-2 

38.95 

15 

1.21E-1 

60.66 

16 

2.92E-1 

145.96 

17 

1.61E-1 

80.56 

7.58E-2 


37.91 


In  order  to  incorporate  a  background  level  into  the 
Bayes'  theorem  analysis,  it  was  necessary  to  generate  a 
separate  background  "source"  based  on  the  high  resolution 
data  used  to  generate  the  original  library  sources  (ref. 
chap.  2).  The  background  measurements  were  taken  at  the  same 
time  the  original  source  measurements  were  made,  and 
therefore  should  accurately  reflect  the  background  in  the  low 
resolution  detector  of  this  problem.  The  relative  intensity 
distribution  of  the  background  is  presented  in  table  V-9.  A 
background  level  of  500  counts  is  assumed,  based  on  the 
background  measurements  taken  with  a  counting  time  of  about 
30  minutes  (as  before).  The  reference  levels  based  on  this 
number  are  also  provided  in  table  V-9.  Since  the  Poisson 
distribution  is  being  used,  the  background  level  will  remain 
constant  at  each  range  (i.e.,  same  counting  time). 

The  reference  level  of  counts  was  added  channel  by 
channel  to  each  reference  source  spectrum,  and  a  random 
(measured)  background  count  level  was  determined  separately 
using  the  same  method  as  was  discussed  in  appendix  G.  Hence 
the  random  number  generator  is  called  upon  twice:  once  to 
generate  the  "normal"  measured  source  spectrum,  due  to  the 
photo  emissions  of  the  nuclear  material,  and  once  to  generate 
the  measured  background  spectrum.  The  results  of  both 
calculations  are  then  combined  to  yield  a  final  measured 
source  spectrum  for  the  analysis  process. 

As  can  be  seen  in  table  V-10  (and  figure  12),  the 
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Table  V-8 


RANGE  Cm 


RANGE  (m 


Bayes'  Posterior  Distribution 
(Multinomial  with  no  background) 
1.  Measured:  Source  A 
P(  A2  I 


50 

1.000 

0.000 

0.000 

100 

1.000 

3. 176E-30 

0.000 

200 

9.999E-1 

8.493E-5 

1 . 101E-11 

300 

2.398E-1 

7.586E-1 

5.949E-4 

400 

3. 743E-1 

4.165E-1 

2.092E-1 

500 

4.452E-1 

3.483E-1 

2.066E-1 

2.  Measured 

Source  B 

RANGE  (m) 

P( A1 | X ) 

PC  A2 1 X ) 

P(A3 I  X) 

2 . 429E-29 

1.470E-1 

1.378E-1 

1.755E-1 

2.673E-1 


1.00 

8.468E-1 

8.622E-1 

3.344E-1 

3.397E-1 


3.  Measured  Source  C 


3 . 852E-14 

6.219E-3 

1.453E-5 

4.902E-1 

3.930E-1 


1 .798E-1 


3.060E-1 


5. 142E-1 


Bayes1  Theorem  With  Multinomial:  No  Background 

As  was  discussed  in  chapter  IV,  the  multinomial 
distribution  is  used  when  the  detector  counts  to  a 
predetermined  number  of  counts,  then  shuts  off.  This  method 
was  programmed  in  order  to  contrast  its  effectiveness  with 
that  of  the  Poisson  distribution,  which  has  no  constraint  on 
the  number  of  counts  collected.  The  posterior  distributions 
from  the  multinomial  runs  are  presented  in  table  V-8. 

Even  though  the  multinomial  analysis  was  done  using  the 
same  initial  reference  count  level  of  1000  counts  at  50 
meters,  it  can  be  seen  that  the  probabilities  of  accurately 
identifying  a  measured  source  drop  off  much  more  drastically 
than  do  the  Poisson  values.  The  assumption  that  each  source 
measured  has  exactly  the  same  number  of  counts  (the  basis  for 
this  distribution)  probably  drives  the  values  down  more 
rapidly  as  the  range  increases.  Recall  in  the  case  of 
Poisson  statistics,  source  A  had  more  total  counts  collected 
at  500  meters  than  did  either  of  the  other  measured  sources, 
which  aided  in  maintaining  a  higher  probability  with 
increased  range. 

The  results  presented  in  table  V-8  are  also  depicted 
graphically  in  figure  11.  Since  the  multinomial  distribution 
yields  lower  values  of  probability,  the  Poisson  methodology 
will  be  used  exclusively  in  the  remaining  analysis  effort. 


Bayes'  theorem  extends  positive  identification  a  minimum  of 
100  meters  for  sources  A  and  C,  and  50  meters  for  source  B 
(reference  table  V.l). 

The  most  striking  aspect  of  the  Bayes'  probabilities  is 
that  they  remain  relatively  high  even  though  the  number  of 
counts  available  is  very  small.  Indeed,  all  three  source 
spectra  look  strikingly  similar  at  400  meters,  yet  measured 
source  A  is  still  identified  78.1%  of  the  time,  and  measured 
source  C  92.6%  of  the  time.  Therefore,  if  the  criterion  for 
positive  identification  is  reduced  to  a  level  less  than  3 
sigma,  Bayes'  theorem  can  extend  the  range  of  identification 
even  further. 

This  is  in  contrast  to  the  method  of  photopeak 
identification,  which  fails  at  relatively  small  ranges 
(around  50  meters).  This  failure  is  due  to  the  small  number 
of  uncollided  photons  available,  since  the  number  is  rapidly 
reduced  by  the  exponential  term  in  Beer's  Law  (appendix  A). 
Even  large  photopeaks  cannot  withstand  the  exponential  losses 
in  conjunction  with  the  rising  background  level  due  to 
downscattered  photons.  Indeed,  in  regions  where  enough 
photons  exist  to  render  photopeak  identification  viable 
(close  to  the  source),  Bayes'  theorem  will  also  return 
favorable  results,  as  evidenced  by  the  100%  values  at  a  range 
of  50  meters. 

The  results  of  the  Poisson  version  of  the  Bayes' 
analysis  are  plotted  in  figure  10. 


the  inability  of  Bayes'  theorem  to  generate  information  on 
non-library  source  composition.  A  possible  solution  to  this 
problem  might  be  to  build  up  a  vast  number  of  library  sources 
based  on  the  gamma  ray  spectra  of  specific  isotopes.  When  an 
unkown  is  encountered,  the  Bayes*  analysis  could  be  run  using 
different  combinations  of  the  isotopic  reference  sources, 
until  a  reasonable  posterior  distribution  is  reached.  The 
probability  of  occurrence,  P(Ai),  could  also  be  modified  in 
order  to  reach  a  high  degree  of  certainty.  Chapter  VI  will 
briefly  discuss  some  follow  on  recommendations  based  on  this 
subject . 


V I .  Summary  and  Recommendations 


To  summarize  this  thesis  problem,  the  purpose  was  to 
test  Bayes'  theorem,  in  theory,  as  an  aid  to  extending  the 
range  of  positive  identification  of  nuclear  materials,  using 
features  in  the  characteristic  gamma  ray  spectra.  To  achieve 
this  goal,  four  basic  steps  were  initiated,  as  covered  in  the 
first  four  chapters  of  this  report.  First,  the  library 
source  spectra  had  to  be  selected,  by  studying  available  high 
resolution  gamma  ray  spectra.  These  were  then  regrouped  into 
a  convenient,  low  resolution  form  for  which  macroscopic  cross 
sections  are  available  (18  groups).  Secondly,  a  transport 
code  was  developed  using  the  Finite  Element  Method  to 
transport  these  library  sources  to  any  desired  range,  in 
order  to  come  up  with  the  reference  sources,  or  in  other 
words,  the  sources  which  might  result  from  the  ideal 
measurement  of  the  source  spectra  at  a  distance.  Third,  the 
reference  sources  were  degraded  into  simulated  measured 
sources  by  applying  various  statistical  distributions,  which 
included  the  Poisson  (and  Gaussian)  and  the  multinomial 
distributions.  And  lastly,  the  measured  sources  were  analyzed 
using  Bayes'  theorem  in  order  to  assess  the  probability  that 
each  measured  source  would  be  correctly  identified.  The 
results  were  compared  with  the  method  of  photopeak 
identification,  to  see  if  any  improvements  could  be  gained  by 
using  Bayes'  theorem. 

As  can  be  seen  by  the  results  in  chapter  V,  Bayes' 


theorem  has  strong  points  and  weak  points.  As  was  seen  in 
the  case  of  measured  sources  with  the  background  extracted, 
Bayes'  theorem  offers  a  dramatic  increase  in  the  accuracy  of 
identification  when  contrasted  with  photopeak  identification. 
When  background  is  not  extracted  from  the  measurements  prior 
to  analysis,  the  results  are  less  convincing  in  the  low  count 
environment  (i.e.,  measurements  taken  at  extreme  distances), 
but  still  offer  an  improvement  over  photopeak  identification. 
The  major  disadvantage  with  Bayes'  theorem  appears  to  be  the 
requirement  for  the  measured  sources  to  be  contained  as 
library  sources  in  order  to  be  identified. 

Recommendations 

Bayes'  theorem  does  indeed  extend  the  range  of  detection 
of  nuclear  materials,  but  much  more  research  needs  to  be 
done.  Some  follow-on  areas  of  research  include  investigating 
the  utility  of  Bayes'  theorem  on  high  resolution  data. 
Theoretically,  the  results  should  be  even  more  dramatic  than 
in  the  low  resolution  case.  This  is  evident  from  equation 
(40)  in  chapter  IV,  where  it  is  seen  that  the  deviations 
between  measured  and  reference  (library)  counts  in  channels 
will  make  the  probability  extremely  small  in  the  case  of  a 
mismatch,  since  the  multiplication  has  increased  from  18 
channels  to  a  much  larger  number  of  channels.  Therefore,  the 
probability  of  measuring  the  wrong  source  will  approach  a 
small  number  (possibly  zero)  much  more  rapidly  than  in  the  18 
channel  scenario.  An  interesting  follow-on  effort  could  test 


Cr 


this  concept  by  increasing  the  number  of  channels  to  some 


(arbitrary)  large  number.  The  18  gamma  group  cross  section 
data  could  be  interpolated  to  accomodate  this  larger  number 
of  energy  groups. 

Secondly,  a  method  should  be  investigated  that  can 
adjust  the  prior  distribution  in  conjunction  with  the  library 
sources  in  order  to  positively  identify  new  sources.  As  was 
mentioned  in  chapter  V,  a  large  number  of  isotopic  library 
sources  might  be  developed  and  used  in  combination  with  the 
prior  distribution  to  determine  if  Bayes'  theorem  can  aid  in 
identifying  new  sources.  As  it  stands  now,  the  results 
discussed  in  chapter  V  show  that  Bayes'  theorem  yields  no 
information  on  new  source  composition,  when  only  a  small 
number  of  library  sources  is  available.  This  aspect  of 
Bayes'  theorem  is  a  potentially  complicated  problem,  and 
unfortunately,  time  didn't  allow  it  to  be  addressed  in  this 
effort . 
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Appendix  A:  Photopeak  Analysis  Considerations 


Photopeak  Analysis  in  Theory. 

As  has  been  mentioned  in  several  places  throughout  the 
main  text,  the  basic  premise  of  this  thesis  is  to  extend  the 
range  of  spectral  identification  by  applying  Bayes'  theorem  to 
measured  spectra  instead  of  the  method  of  photopeak 
identification.  The  principal  difference  in  these  two  methods 
is  that  photopeak  identification  uses  only  the  uncollided 
fluence  as  an  identification  feature,  whereas  Bayes'  theorem 
uses  both  the  uncollided  fluence  and  the  scattered  fluence  as 
useful  information  for  means  of  source  identification. 

Mathematically,  the  virgin  (uncollided)  fluence  is 
attenuated  exponentially,  as  is  seen  in  Beer's  Law: 

I(r)  =  IQ(r)  exp(-Ztr)/4irr  2  (44) 

In  equation  (44),  Ztis  the  total  macroscopic  cross  section  for 
interaction  for  photons  at  the  given  energy.  Hence,  any 
interaction,  including  absorption  and  scattering,  removes  a 
photon  from  the  uncollided  fluence  catagory. 

In  order  for  a  photopeak  to  be  useful  for  source 
identification  purposes,  it  must  be  a  minimum  of  three  standard 
deviations  (sigma)  above  the  background  level.  In  the  case  of 
the  Poisson  distribution,  one  standard  deviation  is  the  square 
root  of  the  number  of  counts  in  the  given  channel.  A  three 
sigma  level  provides  a  99.73  percent  certainty  of  positive 


identification  (4).  Therefore,  when  a  photopeak  falls  below 
the  three  sigma  level,  the  method  of  photopeak  identification 
is  no  longer  reliable.  Since  Bayes’  theorem  uses  a  portion  of 
the  downscattered  photons  for  identification  purposes,  it 
should  in  theory  extend  the  range  of  source  identification 
beyond  the  range  where  the  photopeaks  drop  below  the  three 
sigma  level. 

Photopeak  Analysis  in  Practice. 

In  this  problem,  the  three  sigma  criterion  must  be 
evaluated  at  each  range  at  which  the  Bayes'  theorem  analysis  is 
being  performed.  This  involves  first  selecting  distinct 
photopeaks  in  each  high  resolution  source  which  makes  each 
source  unique. 

Each  of  the  original  high  resolution  sources  was  examined 
for  one  very  distinct  line  for  the  purpose  of  comparison.  In 
source  A  the  line  occurs  in  channel  10,  in  source  B  it  occurs 
in  channel  12,  and  in  source  C  the  "test"  photopeak  also 
appears  in  channel  12.  These  photopeaks  were  then  degraded 
using  the  following  assumptions: 

(1)  the  line  occurs  at  the  center  of  the  low 
resolution  channel; 

(2)  the  line  is  infinitesimally  narrow,  so  that 
no  photons  can  downscatter  into  the  photopeak; 

(3)  the  background  (continuum)  in  the  low 
resolution  bin  is  divided  equally  among  all 

the  high  resolution  channels  present  in  the  bin. 
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In  order  to  examine  the  "best  case"  for  the  method  of 
photopeak  identification,  the  relative  intensity  of  each  of 
the  photopeaks  at  the  source  is  used  to  calculate  the  number 
of  uncollided  counts  present  in  the  photopeak  at  a  range  of 
50  meters.  The  count  levels  in  the  low  resolution  channels 
are  obtained  from  tables  V.2  through  V.7  (Poisson 
Distribution  with  no  background).  The  number  of  counts  in 
each  photopeak  is  determined  using  Beer's  Law.  This  number 
is  then  subtracted  from  the  number  of  measured  counts  in  the 
appropriate  low  resolution  channel,  and  the  rest  of  the 
counts  are  divided  equally  among  the  high  resolution 
channels,  giving  a  smooth  background  level  in  accordance  with 
assumption  (3)  above.  Figure  13  illustrates  the  relationship 
between  the  low  resolution  bin  and  the  high  resolution 
structure . 

As  can  be  seen  in  figure  13,  as  the  range  increases,  the 
photopeak  will  degrade  in  accordance  with  Beer's  Law  and  the 
relative  background  level  will  increase  as  photons 
downscatter  into  the  bin  from  higher  energy  channels. 

The  analysis  was  performed  on  the  three  photopeaks 
discussed  above,  and  the  results  are  contained  in  table  A-l. 


Number  of  Counts 


TABLE  A- 1 


Photopeak  Degradation 

1.  Source  A 

COUNTS 

50  meters 

100  meters 

200  meters 

ft 

total: 

6.00 

1.00 

1.00 

■*-  V 

photopeak : 

4.00 

8.27E-9 

3.87E-10 

y;/ 

« y 

sigma : 

2.39E-1 

1.69E-1 

1.69E-1 

1 

3  sigma: 

7.17E-1 

5.07E-1 

5.07E-1 

l‘X" 

bkgd .  : 

5.71E-2 

2. 

2.86E-2 

Source  B 

2.86E-2 

L 

fv< 

COUNTS 

50  meters 

100  meters 

200  meters 

total : 

7.00 

2.00 

0 

photopeak : 

1.30 

2.38E-9 

0 

ft 

sigma : 

5.34E-1 

3.16E-1 

0 

*  *  V 

3  sigma: 

1.60 

9.49E-1 

0 

y’s-i 

V/V 
*  *  “J 

bkgd . : 

2.85E-1 

3^_ 

1 .00E-1 

Source  C 

0 

m 

ITT 

ftY* 

,'V*- 

COUNTS 

50  meters 

100  meters 

200  meters 

total : 

1 

0 

0 

ft 

photopeak : 

1 

0 

0 

W* 

sigma : 

0 

0 

0 

3  sigma: 

0 

0 

0 

ft 

bkgd  : 

0 

0 

81 

0 

1 

i 

1 

x 
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Appendix  B:  18  Gamma  Group  Cross  Sections 


This  appendix  contains  the  data  file  listing  of  the  18 
group  sea  level  air  gamma  cross  sections  used  in  the 
radiation  transport  process.  All  cross  sections  are 
macroscopic,  the  units  being  inverse  centimeters.  The  list 
begins  with  gamma  group  one  (highest  energy  group)  and 
continues  down  to  group  eighteen.  The  first  two  numbers  in 
each  group  (separated  by  a  comma)  represent  the  group 
transport  cross  section  and  the  group  total  cross  section,  i 
that  order.  The  next  value  in  each  group  is  the  cross 
section  for  scatter  from  group  one  into  the  current  group. 
The  next  cross  section  is  the  scatter  from  group  two,  and  so 
on,  until  within  group  scatter,  which  is  the  last  cross 
section  in  each  group  (14). 


2.6781E-5.3.0663E-5 

1.2991E-6 

2.9082E-5.3.3578E-5 
2. 1316E-6 
1 .5056E-6 

3.0430E-5.3.7458E-5 
2 . 2584E-6 
3 . 2793E-6 
2 . 3603E-6 

3.4609E-5.4.2309E-5 
1 .6760E-6 


(Group  1,  8000  -  10000  KeV) 


(Group  2,  6500  -  8000  KeV) 


(Group  3,  5000  -  6500  KeV) 


(Group  4,  4000  -  5000  KeV) 
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3.4861E-6 
2 . 5878E-6 

3.5861E-5.4.8296E-5 
1 .9585E-6 
2 . 6055E-6 
3 . 7052E-6 
5 . 6385E-6 
4.2018E-6 

4.4473E-5.5.4819E-5 
1 .  1661E-6 
1 . 5069E-6 
2.0530E-6 
2 . 9597E-6 
4 . 6506E-6 
3.4867E-6 

4.5947E-5.6. 1097E-5 
1.3771E-6 
1 . 7477E-6 
2.31 40E-6 
3.2071 E-6 
4.7851 E-6 
7 . 3590E-6 
5. 1337E-6 

5.2539E-5.6.8224E-5 
1  .  1220E-6 
1 .4045E-6 


(Group  5,  3000  -  4000  KeV) 


(Group  6,  2500  -  3000  KeV) 


(Group  7,  2000  -  2500  KeV) 


(Group  8,  1660  -  2000  KeV) 


N 


N 


a 
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1 .8185E-6 


2 . 4386E-6 
3 . 4698E-6 
5 . 0735E-6 
7 . 4414E-6 
5.31 79E-6 

5.3610E-5, 7.5903E-5  (Group  9,  1330  - 

1 . 3140E-6 

1 .6295E-6 

2 . 0759E-6 

2 . 7 143E-6 

3.7127E-6 

5. 1831E-6 

7 . 2570E-6 

1 . 0756E-5 

7.6158E-6 

5. 1254E-5.8.6398E-5  (Group  10,  1000 

1 .6733E-6 

2.0597E-6 

2 . 5888E-6 

3 . 3107E-6 

4.3618E-6 

5 . 8008E-6 

7 . 6914E-6 

1 .0725E-5 

1 .5861E-5 


1660  KeV) 


-  1330  KeV) 


1 .2197E-5 


6. 1939E-5.9.8011E-5 
1 . 8052E-6 


(Group  11,  800  -  1000  KeV) 


1 .6008E-6 
1 .9948E-6 
2.5154E-6 
3 . 2290E-6 


4. 1378E-6 


5 . 2363E-6 
6 . 8793E-6 


9 . 4990E-6 


1.5564E-5 


1 .2529E-5 


5.4339E-5, 1 . 1023 E-4 
1 . 6878E-6 


(Group  12,  600  -  800  KeV) 


2 . 0662E-6 
2 . 5674E-6 


3.2143E-6 


4 . 0659E-6 


5 . 0889E-6 


6 . 2289E-6 


7. 7982E-6 


1 . 0093E-5 


1 .5022E-5 


2.5010E-5 


1 . 9937E-5 


3.2786E-5, 1 .2779E-4 


(Group  13,  400  -  600  KeV) 


1 .5137E-5 


1 .3442E-5 


1 . 1976E-5 
1 .0767E-5 
9 . 9962E-6 
9. 7302E-6 
1 .0069E-5 
1 . 1029E-5 
1 . 2672E-5 
1 .6497E-5 
2.3561 E-5 
3.8612E-5 
3 . 6 189E-5 

4.9010E-5, 1 .4714E-4  (Group  14,  300  -  400  KeV) 

1 . 7477E-6 

2. 1519E-6 

2 . 6905E-6 

3.3851 E-6 

4. 2817E-6 

5 . 308 1 E-6 

6. 3437E-6 

7. 5724E-6 

8 . 986  1  E-6 

1 . 1101E-5 

1 .3958E-5 

1 .8611E-5 

3 . 538 1 E-5 


3. 7459E-5 


2/2 


AD-A159  244 
UNCLASSIFIED 


IDENTIFICATION  OF  NUCLEAR  MATERIALS  FROM  REMOTE 
DETECTION  OF  CHARACTERIST  <U)  AIR  FO*CE  INST  OF  TECH 
URIGHT-PATTERSON  AFB  OH  SCHOOL  OF  ENG I  L  U  BRASURE 
MAR  85  AFIT/GNE/PH/85M-2  F/G  18/4  NL 


MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BUREAU  OF  STANDARDS-1963-A 


1 . 7734E-5 , 1 . 6655E-4 
1 . 1719E-6 


(Group  15,  200  -  300  KeV) 


1 .5017E-6 
1 . 9791E-6 
2 . 6574E-6 
3  *  65 1 5E-6 
4.9480E-6 
6.4615E-6 
8.5184E-6 
1 . 1264E-5 
1 . 6018E-5 
2 . 2300E-5 
2 . 7297E-5 
3 . 6680E-5 
6 . 4487E-5 
6 . 4947E-5 


5.2605E-5.1. 9710E-4 


(Group  16,  100  -  200  KeV) 


6 . 5306E-7 


5 . 7746E-6 
1 . 9534E-5 
4 . 5196E-5 
1 . 0160E-4 
1 .4823E-4 

1.4686E-4,2. 4089E-4 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 

4 . 8024E-5 
2.0558E-4 

3.6340E-4.4.0987E-4 


(Group  17,  50  -  100  KeV) 


(Group  18,  20  -  50  KeV) 


Appendix  C:  Simplifying  the  Three-Dimensional 
Transport  Problem 


This  appendix  will  outline  the  steps  by  which  the 
three-dimensional,  time  independent,  radially  symmetric 
diffusion  equation  is  reduced  into  a  one-dimensional  problem. 
By  simplifying  the  diffusion  equation  in  this  manner,  the 
approximating  polynomials  developed  for  the  one-dimensional 
problem  in  chapter  III  are  still  applicable. 

Beginning  with  the  time  independent,  spherically 
symmetric  form  in  three  dimensions 


i* 


-D  7-  -h  ^  -  s  (45) 

Next,  a  simple  substitution  is  made  in  equation  (45),  as 
follows 


u ( r )  =  rF( r )  (46) 

By  solving  equation  (46)  for  F,  substituting  back  into 
equation  (45)  and  collecting  terms,  the  resulting  equation  is 

_  D  iiu  +  £a—  =  S  (47) 

r  dr  2  r 
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Which  simplifies  to  the  form 


_D  lis.  +  zau  =  s,  (48) 

dr'' 

where  S  '  =  rS  (49) 

Equation  (48)  is  in  a  form  similar  to  the  one  dimensional 
diffusion  equation  (chapter  III,  equation  (26))  and  can  be 
solved  using  the  same  basic  code,  with  the  modification  to 
the  source  term  given  by  equation  (49).  Hence  each  group  in 
the  multigroup  diffusion  problem  is  solved  in  terms  of  the 
variable  u,  and  the  fluence  at  each  node  is  obtained  by 
applying  equation  (46). 


Appendix  D:  The  Analytical  Solution 
(Simple  Case) 

In  order  to  verify  the  radiation  transport  code,  it  is 
necessary  to  look  at  a  simple  case  of  the  diffusion  equation 
and  obtain  a  solution  analytically.  Hence,  a  simple  two 
energy  group  problem  will  be  considered,  and  the  diffuse 
solution  for  the  lower  energy  group  will  be  solved. 

First,  the  homogeneous  solution  to  equation  (48)  in 
appendix  C  must  be  found 

-D  +  I  u  =  0  (50) 

dr  2 

The  homogeneous  solution  follows  quite  easily,  and  by 
requiring  the  solution  to  remain  finite  with  increasing 
range,  r,  it  reduces  to  one  term. 

u  (r)  =  C  exp( -( Ea/ D ) ar )  (51) 

h 

Next,  the  particular  solution  to  equation  (48)  is  determined 
by  first  assigning  u  as 

u  (r)  =  A  r  +  A  (52) 

p  1  U 

By  differentiating  equation  (52)  twice  with  respect  to  r  and 
substituting  back  into  equation  (48)  from  appendix  C,  the 


'wn  wv  1 


r 
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r 

K 
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result  becomes 


IaA1r  +  ZaA  Q  =  rS  (53> 

By  equating  the  like  terms  of  equation  (53),  the  coefficients 
are  determined  as  follows 


A  0  =  0  (54) 

kY  =  S/Za  (55) 


The  final  solution  becomes 

u(r)  -  C  exp(-(Za/D)  *r)  +  S/Za  r  (56) 

To  determine  the  constant  C,  the  first  boundary  must  be 
applied.  Since  there  is  only  virgin  flux  at  the  origin  of 
the  source,  the  diffuse  solution  must  be  zero  when  r  is  zero, 
and  hence 


u  ( 0 )  =  C  =  0 


(57) 


Equation  (56)  now  reduces  to  one  term,  and  when  it  is  solved 
for  the  fluence  by  applying  equation  (46)  (appendix  C),  it 
becomes 


i* 


F(r)  =  S(r)/Za 


(58) 
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The  initial  intensity  levels  assigned  to  these  two 
energy  groups  were  0.40  counts  per  photon  (group  1)  and  0.35 
counts  per  photon  (group  2).  The  values  of  the  diffusion 
constant  and  macroscopic  absorption  cross  section  were 
2.6637E+4  cm,  5 . 8079E-5/cm ,  2.6232E+4  cm  and  5 . 5647E-5/cm ,  in 
that  order.  The  results  for  the  diffuse  solution  using  the 
analytical  solution  was  4.01E-11  counts  per  photon  per  square 
cm,  and  the  transport  code  solution  for  only  one  mesh  space 
was  8.65E-12  counts  per  photon  per  square  cm. 


Appendix  E:  Transport  Code  Listing 


> 


This  appendix  provides  a  listing  of  the  radiation 
transport  code  discussed  in  chapter  III.  The  subroutines  are 
discussed  one  at  a  time,  along  with  a  brief  description  of 
the  variables.  Variables  which  are  discussed  in  the  main 
program  will  not  be  discussed  a  second  time  in  each 
subroutine . 

When  running  this  code,  it  is  very  important  to  make 
dimensional  changes  to  two  of  the  variables,  in  order  to 
maintain  the  positive  definite  nature  of  the  arrays  which  are 
processed  by  the  IMSL  library  subroutines.  These  include  the 
S  array  in  the  main  program  and  the  S  and  the  MG  arrays  in 
subroutine  diff.  The  S  array  must  be  dimensioned  as  a  NV  * 
3  matrix  and  MG  as  a  NV  *  4  matrix.  Here  NV  stands  for  the 
number  of  variables,  and  is  equal  to  two  times  the  number  of 
nodes.  As  listed,  the  program  is  set  up  for  100  mesh  spaces, 
or  202  variables  (the  flux  and  current  at  each  node).  If  the 
number  of  mesh  spaces  is  changed  without  redimensioning  these 
two  arrays  accordingly,  the  IMSL  subroutines  will  return 
error  messages. 

Another  important  note  concerning  this  program  concerns 
the  lines  containing  a  "D"  where  a  comment  code  "C"  would 
normally  go.  These  lines  are  debugging  lines,  and  are 
compiled  on  the  Harris  800  when  the  debug  mode  is  used, 
otherwise  they  are  treated  as  comment  lines  (not  compiled). 

This  code  is  written  in  Fortran  77,  and  extensive  use 
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was  made  of  references  (1)  and  (9),  as  well  as  Dr.  Donn 
Shankland's  knowledge  of  the  operating  system.  The  IMSL 
documentation  was  also  extremely  helpful  for  solving  the 
matrix  equations  in  this  program. 

Program  Description 
Main  Program. 

The  main  program  is  responsible  for  calling  all  the 
subroutines  and  finally  printing  out  the  final  solution  for 
the  fluence  at  each  node. 

List  of  Variables. 

R  -  stores  the  ranges  from  the  source, 

XTR  -  stores  the  transport  cross  sections, 

XR  -  stores  the  removal  cross  sections, 

XT  -  stores  the  total  cross  sections, 

XS  -  stores  the  scattering  cross  sections, 

SG  -  stores  one  group's  removal  cross  sections, 

C  -  stores  the  boundary  conditions  cl,c2,dl,d2, 

B  -  stores  the  boundary  conditions  el,e2 

DD  -  stores  the  diffusion  constants  for  one  group, 

S  -  stores  the  source  matrix  from  Simpson's  approx., 
C1,C2,D1,D2  -  boundary  conditions, 

E1,E2  -  boundary  conditions, 

T  -  stores  values  of  1,  2*PI  &  4*PI, 

MAXN  -  maximum  number  of  nodes, 

MAXMAT  -  maximum  number  of  materials  (1  in  this  prob.) 


N  -  number  of  nodes. 


K  -  geometry  factor, 

MAT  -  stores  the  material  number, 

NMAT  -  number  of  materials, 

NC  -  IMSL  subroutine  variable, 

I  -  integer  counter, 

NV  -  number  of  variables, 

MAXG  -  number  of  groups  used. 

Subroutine  GDATA. 

This  subroutine  retrieves  the  basic  data  for  setting  up 
the  appropriate  mesh  spacing. 

List  of  Variables. 

A,  AA,  B,  BB,  C  -  position  varibles, 

NR  -  number  of  regions, 

MN  -  material  number, 

NS  -  number  of  mesh  spaces  per  region, 

M  -  position  variable. 

Subroutine  MDATA. 

This  subroutine  retrieves  data  on  each  different  group 
and  each  different  material. 

List  of  Variables. 

I,  J,  K  -  integer  counters. 

Subroutine  BDATA, 

This  subroutine  retrieves  the  boundary  data  and  sets  up 
the  C  matrix. 

Subroutine  BOUNDR, 

This  subroutine  computes  the  boundary  condition  matrix, 
C,  in  the  case  of  the  radiation  transport  problem,  due  to  the 


energy  dependent  boundary  conditions. 


Subroutine  DIRECT. 

This  subroutine  calculates  the  direct  fluence  at  the 
nodal  points  and  at  the  midpoints  of  each  mesh  interval, 
values  which  are  modified  in  subroutine  TERM  and  ultimately 
used  in  the  Simpson's  approximation. 

List  of  Variables. 

GROUP  -  the  current  group  number, 

IA  -  the  current  interval  number  (mesh  space), 

GAMMA  -  stores  the  initial  intensity  levels  of 
all  the  groups, 

DPHI  -  stores  the  direct  fluence  at  each  node, 

MPHI  -  stores  the  direct  fluence  at  each 
half  node, 

MULT  -  temporary  storage  variable, 

MIDPT  -  the  midpoint  of  the  interval. 

Subroutine  TERM. 

This  subroutine  modifies  the  direct  fluences  into  the 
source  term  used  by  the  Simpson's  approximation. 

List  of  Variables. 

IH  -  group  counter  for  downscatter  calculations, 

MSOURCE  -  source  terra  at  mid-interval, 

DR  -  distance  from  the  left  node  to  the 
midpoint  of  the  current  interval. 

Subroutine  DIFF . 

This  subroutine  sets  up  and  solves  the  matrix  equation 
derived  in  chapter  III.  In  arriving  at  the  solution  for  the 


salve  far  lambdas 


C 

c 
c 

0  PRINT*," 

D  PRINT*, 'entering  LEQT1P  to  solve  for  lamdas' 

CALL  LEQT1P  (Q,1 ,2,R,2,ID,D1 ,02, IR) 

D  PRINT*,'  solution  follows  (lamdas):' 

D  PRINT*,'  LAMOA  1  =' ,R(1 ) 

D  PRINT*,'  LAND A  2  =',R(2) 

C 

C  add  boundary  terms  to  the  free  solution 

C 

D  PRINT*,'  entering  SAXPY  to  add  in  boundary  terms  . 
DO  7  I  =  1 ,2 

CALL  SAXPY  (NV ,R( I )  ,S(1 ,1+1 ),1,S,1) 

7  CONTINUE 

D  PRINT*," 

D  PRINT*,'  the  solution  for  the  variable  U  follows' 

D  PRINT*,'  (for  fluence*r  AND  current*r)' 

D  PRINT*,'  (still  in  subroutine  DIFF)' 

D  DO  500  IP  =  1  ,NV 

D  PRINT*, S(IP,1  ),S(IP,2),S(IP,3) 

D500  CONTINUE 

C  since  the  solution  is  in  terms  of  the  arbitrary 

C  variable  "U",  which  is  equal  to  R  *  FLUENCE, 

C  we  must  now  SOLVE  far  the  scattered  fluence  by 

C  dividing  "U"  by  R  at  each  node  and  putting  the 

C  solution  back  in  the  matrix  "S". 

DO  600  I  =  1 ,NV,2 

IX  =  I  -  (1/2  -  0.5) 

XX(1)  =  1 

D  PRINT*, 'S( ' ,1, ' ,1 )  =' ,S(I,1 ) 

D  PRINT*, 'XX( ' , IX,')  =',XX(IX) 

3(1,1 )  =  S(I ,1 )/XX(IX) 

D  PRINT*, S(1, 1  ) 

600  CONTINUE 


soli/e  the  main  matrix  problem 


C 

c 
c 

D  PRINT*," 

D  PRINT*,'  entering  LEQ1P8  .  . 

D  PRINT*," 

CALL  LEQ1PB  (l*IG,NV,3,NV,S,NV,3,IDGT,D1  ,D2,IER) 

0  PRINT*," 

D  PRINT*,'  the  solution  for  tO,  tl  &  t2  follows 

0  00  138  IY  =  1 ,NV 

0  PRINT*,S(IY,1 ),S(IY,2),S(IY,3) 

D138  CONTINUE 

C 

C  compute  the  small  matrix 

C 

K=0 

0  PRINT*," 

0  PRINT*, 'entering  VIPRFF  imsl  subroutine  now...' 

D  PRINT*," 

D  PRINT*," 

D  PRINT*, 'the  pre-VIPRFF  C  matrix  follows:' 

0  00  145  IT  =  1 ,NV 

0  PRINT*, C(IT,1 ),C(IT,2) 

D1A5  CONTINUE 

DO  200  I  =  1 ,2 

CALL  VIPRFF  (C(1,I),S,NU,1,1,R(I)) 

0  PRINT*," 

0  PRINT*, 'for  i  =',I 

0  PRINT*, 'R(', I,')  =' ,R(I) 

0  PRINT*, 'B(', I,')  = '  ,8( I ) 

R(I)  =  8(1)  -  R( I ) 

0  PRINT*, 'B(I)  -  R( I )  =' ,R(I) 

DO  150  J  =  1,1 

K  =  K  ♦  1 

D  PRINT*," 

D  PRINT*,'  entering  VIPRFF  again  ...' 

CALL  VIPRFF(C(1 ,1) ,S(1 ,J+1 ),NU,1,1 ,Q(K)) 

D  PRINT*,'  for  j=',J 

0  PRINT*,'  and  k  =',K 

0  PRINT*,'  Q(',K,')  =' , 0 ( K ) 

150  CONTINUE 
200  CONTINUE 


c 

C  compute  the  source  term 

C 

C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  add  in  the  downscatter  contribution  ...  C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

DO  80  IH  =  1,  IG-1 

CALL  TERM ( N .MAT , IG , IH , XX , XS , PHI , SOURCE , 

+  MPHI.MIDPT.IA.MSOURC.XT) 

C  return  with  "SOURCE"  array  values  and  compute 

C  the  source  matrix  using  Simpson's  Approximation. 

C  NOTE:  downscattered  sources  are  added  in  to  the 

C  inscattered  source  already  present  in  the  "S"  array. 

SA  =  H  *  S0URCE(IA)/6 
SB  =  H  *  P1S0URC  /  3 
SC  =  H  *  S0URCE(IA+1 )/6 
S(K+1 ,1 )  =  S(K+1 ,1 )  +  SA  +  SB 
S(K+2,1)  =  S(K+2,1)  -  .25*H*SB/D 
S(K+3,1)  =  S(K+3,1)  +  SB  +  SC 
S(K+4,1)  =  S(K+4,1)  +  .25*H*SB/D 
80  CONTINUE 

C 

C  insert  boundary  conditions 

C 

DO  90  J=1 ,2 

S(K+1 ,  J+1 )  =  C(K+1,J) 

S(K+2,J+1)  =  C(K+2, J) 

90  CONTINUE 

100  CONTINUE 
IA  =  NI  +  1 
K  =  IA  +  IA  -  2 
00  125  J  =  1,2 

S(K+1,J+1)  =  C(K+1,J) 

S(K+2,J+1 )  =  C(K+2, J) 

125  CONTINUE 
0  PRINT*," 

0  PRINT*,'  the  MG  matrix  follows  (banded  form)  .  .  .' 

0  DO  135  IZ  =  1,NV 

0  PRINT*, MG(IZ,1),MG(IZ,2),MG(IZ,3),MG(IZ,4) 

0135  CONTINUE 
D  PRINT*," 

0  PRINT*, 'the  source  matrix  follows  (s,c1,c2)  .  .  .' 

D  DO  136  IZ  =  1 ,NV 

D  PRINT *, S ( IZ,1 ),S(IZ,2),S(IZ,3) 

01 36  CONTINUE 


no 


c 

c 

c 


now,  build  each  local  matrix  (for  each  interval) 


(MI  =  N-1 

D  PRINT*,  'N  =',N 

D  PRINT*, 'NV  =  ',NV 

D  PRINT*, 'NI  =  ',NI 

DO  100  IA  =  1 ,NI 
D=DD(IA) 

H=XX(IA+1 )  -  XX(IA) 

D  PRINT*,' for  interval  number', IA 

D  PRINT*,'  D  =  ' ,D 

0  PRINT*,'  XX(IA+1 )  =' ,XX(IA+1 ) 

D  PRINT*,'  XX(IA)  =  ' ,XX(IA) 

D  PRINT*,'  H  =  '  ,H 

D  PRINT*, 'SG(IA)  =' ,SG(IA) 

D  PRINT*," 

r(L(1,1,IA)=1.2*D/H  +  13.*SG(IA)*H/35. 

NL(2,1,IA)=-0.1  -  11  ,*SG(IA)*H*H/(210.*D) 
NL(3,1,IA)=-1.2*0/H  +  9.*SG(IA)*H/70. 
m.(4,1 ,IA)=-0.1  +  13.*SG(IA)*H*H/(420.*D) 
nL(2,2,IA)=(D/(30.*H)+H*SG(IA)/420.)  *  (2.*H/D)**2 
PTU  (3,2,IA)  =— WL  ( 4,1  ,IA) 

NL(4,2,IA)=-(D/(30.*H)+H*SG(IA)/140.)  *  (H/D)**2 
™.(3,3,IA)=rnL(1,1,IA) 

ML(4t3,IA)=-ML(2,1,IA) 

rnL(4,4,IA)=WL(2,2,IA) 

C  now  fill  in  the  rest  of  the  elements 

DO  30  I  =  1,3 
IP  =  I  +  1 
DO  20  J  =  IP, 4 

I¥1L(I,J,IA)  =  NL(J,I,IA) 

20  CONTINUE 

30  CONTINUE 

0  PRINT*," 

D  PRINT*,'  echo  the  local  matrix  now  .  .  .' 

D  DO  35  IZ  =  1,4 

D  PRINT*, 1*L(IZ,1  ,IA),TO.(IZ,2,IA),nL(IZ,3,IA),(»L(IZ,4,IA) 

035  CONTINUE 

C 

C  assemble  the  global  matrix  now 

C 

K=IA+IA-2 
DO  50  I  =  1,4 
DO  40  J  =  1,1 

MG(K+I,4+J-I)=NG(K+I,4+J-I)+I*L(I,J,IA) 

40  CONTINUE 

50  CONTINUE 

D  PRINT*," 

D  PRINT*,' - 

D  PRINT*,'  GROUP  \IG,'  INTERVAL',  I A 

0  PRINT*,' - 

0  PRINT*," 
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*  SUBROUTINE  OIFF  * 

»»##»»***# a*##*#**#*##**#**##**#**#*####*#***#***#*#*##**#***##*###*#*#**# 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

D 

D 

D 

0 

DIO 

C 


C 

C 


4 

C 


3 


SUBROUTINE  DIFF(S,C,B,XX,DD,SG,NV,N,NG,XS,XT,MAT,IG,T) 

INTEGER  NV,N,NI,IG,IH, INSCATTER, NG,MAT(*),K,IA 
REAL  0,DD(*) ,H,SG(*) ,SA, SB, SC, SOURCE (101 ) »T(3) 

+  ,D1,D2,MG(22,4),ML(4,4,100),S(22,3),C(202,2) 

REAL  B(*) ,R(2) ,Q(3) ,XX(101 ) ,XS(5,21 ,21 ) ,DPHI(101 ) , XT (5,21 ) 
+  ,MPHI,MIDPT,MS0URC,PHI(18,101 ) 


This  subroutine  uses  the  FINITE  ELEMENTS  METHOD  to 
solve  the  one-dimensional  diffusion  equation. 


MA  =  1 / ( 30*H )  * 

3B 

-3H/D 

(2H/D)**2 

-36 

3H/D 

36 

-3H/D 

-(H/D)**2 

3H/D 

(2H/D)**2 

MB  =  H/420  * 

15B 

-22H/D 

(2H/D)**2 

54 

-13H/D 

156 

13H/D 

-3(H/D)**2 

22H/D 

(2H/D)**2 

note:  these  matrices  come  from  the  quadratic  terms  of  the 
penalty  function  and  are  symmetric  (hence  the  blank 
entrees  in  the  upper  portions). 

m  =  2  »  n 

PRINT* , ' echo  of  C  matrix  within  subroutine  diff...' 

PRINT*," 

DO  10  I  =  1,NV 

PRINT*, C(I,1),C(I,2) 

CONTINUE 

ZERO  THE  ARRAYS 
DO  3  I  =  1,NV 
DO  4  3  =  1,3 

ZERO  OUT  THE  GOBAL  MATRIX 
MG(I,J)  =  0.0 

ZERO  OUT  THE  SOURCE  MATRIX 
S(I,J)  =  0.0 
CONTINUE 

ZERO  OUT  THE  LAST  COLUMN  OF  THE  GLOBAL  MATRIX 
MG( 1,4)  =  0.0 
CONTINUE 
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*  SUBROUTINE  TERM  * 

it*#*#**##*###*##*****##**#***#****#**##**#*#*#**#*##)*#********#******#**** 


SUBROUTINE  TERM  (N, MAT, IG,IH,R,XS,PHI, SOURCE, 

+  MPHI , MIDPT, I A, MSOURC,  XT ) 

REAL  R( 1 01 ),XS{5,21 ,21 ) ,PHI(1 8,101 ) ,SOURCE(») ,MPHI,MIDPT, 

+  MS0URC,DR,XT(5,21 ) 

INTEGER  N,MAT(*),IG,IH,IA 

C  this  subroutine  computes  the  source  term  which  is  used 

C  in  the  Simpson's  Rule  Approximation  for  the  integrand 

C  term:  -2FSprime. 

C  NOTE:  since  Sprime  =  rS,  the  source  term  becomes: 

C  XS  *  R  *  DPHI 

C 

D  PRINT*," 

D  PRINT*,' in  subroutine  TERM 

0  PRINT*," 

0  PRINT*, 'IG  =',IG,'IH  =' ,IH 

DO  200  I  =  1.N-1 

SOURCE ( I )  =  XS(MAT(I),IG,IH)*PHI(IH,I)*R(I) 

200  CONTINUE 

SOURCE(N)  =  XS(MAT(N-1),IG,IH)*PHI(IH,N)*R(N) 

D  PRINT*, '(IG  STILL  \IG,')  IH  =',IH 

D  PRINT*, 'XS(IG.IH)  =' ,XS(MAT(N-1 ) ,IG,IH) 

C  now  calculate  the  4(pi)r**2  fluence  at  the  left  side  of 

C  this  interval  ... 

MIDPT  =  0.5  *  (R(IA)  +  R(IA+1 ) ) 

MPHI  =  A  *  3.141 5926  *  (R(IA)**2)  *  PHI(IH,IA) 

DR  =  MIDPT  -  R(IA) 

0  PRINT*, 'R(IA)  =' ,R(IA) , '  MIDPT  =',MIDPT,'  DR  =',DR 
D  PRINT*,'  IA  =’ ,IA, 1  4  IA+1  =',IA+1 

C  calculate  the  fluence  at  mid-interval  using 

C  exponential  loss  by  absorption  ONLY  (sigma  removal). 

MPHI  =  MPHI  *  EXP (-XT (MAT( IA) ,IH)*0R)  /  (4*3.1 41 593*MIDPT**2 ) 
MSOURC  =  XS(MAT(IA) ,IG ,IH)  *  MPHI  *  MIDPT 

D  PRINT*,'  MSOURC  =', MSOURC 

D  PRINT*,'  SOURCE(IA)  =' ,S0URCE(IA) 

D  PRINT*,'  SOURCE (IA+1 )  =' ,S0URCE(IA+1 ) 


assign  the  geometry  factor  "K"  in  next  line. 


imJLT  =  GAFFIA(GROUP) 

iterate  over  all  the  nodes 
DO  100  I  =  1 ,N 

IF  (I.EQ.1)  THEN 
R(I)  =  1 

i.e.,  set  the  left  boundary  =  1  cm 
DPHI ( 1 )  =  F!ULT/(T(K)*R(1)**K) 

ELSE 

H(I)  =  R(I)  -  R(I-1) 

MULT  =  MULT*EXP(-XT(FIAT(I-1),GR0UP)»H(I)) 
OPHI(I)  =  FIULT/(T(K)*R(I)**K) 

ENDIF 

CONTINUE 

ends  loop  over  all  nodes 


now  calculate  the  third  direct  fluence  value,  that  is, 
the  value  in  the  middle  of  the  interval  of  interest. 


FIIOPT  =  0.5  *  (R(IA)  +  R(IA+1 ) ) 

FIPHI  =  GAMMA (GROUP)  *  EXP ( -XT ( MAT ( I A ) ,  GROUP )  *1*11  DPT ) 
/  (T(K)  *  MIDPT**K) 

PRINT*,  'DPHI  at  IA  =\IA,'  is  =',0PHI(IA) 

PRINT*, 'DPHI  at  IA+1  =',IA+1,'  is  =',DPHI(IA+1 ) 
PRINT*,  'MPHI  ='  ,MPHI 
PRINT*,'MIDPOINT  =' ,MI0PT 
END 
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*  SUBROUTINE  DIRECT  * 

»*#»*####****#»*»*#*#*»**«#***#*#*###**#***»****#***#*##*******#*#*#**# 

SUBROUTINE  DIRECT  (N,T,R, WAT, XT, DPHI, GROUP, IA,MPHI,MIDPT) 

C 

C  this  subroutine  computes  the  direct  fluence  at  each 

C  node  for  each  group,  filling  the  "DPHI"  array. 

C 

REAL  R(*),XT(5,21),DPHI(101),GAMMA(18),T(3),MULT,H(101  ) 

+  .mPHI.miDPT 
INTEGER  N,K,MAT(»), GROUP, IA 
0  PRINT*," 

D  PRINT*,'  in  subroutine  DIRECT  .  . 

D  PRINT*," 

D  PRINT*, 'GROUP  =', GROUP 

GAMMA (1 )  =  0.0 
GAMMA(2)  =  0.0 
GAMMA(3)  =  0.0 
GAMMA(4)  =  0.0 
GAMMA ( 5 )  =  0 
GAMMA(6)  =  0.33 
GAMMA (7)  =  0.13 
GAMMA (0)  =  0.05 
GAMMA (9)  =  0.10 
GAMMA(IO)  =  0.25 
GAMMA(II)  =  0.01 
GAMMA(12)  =  0.03 
GAMMA(13)  =  0.03 
GAMMA(U)  =  0.04 
GAMMA(15)  =  0 
GAMMA(16)  =  0.03 
GAMMA(17)  =  0 
GAMMA(IB)  =  0 


***•**•«•*************«*«*****»••***»*«*******•**»•*****•»»«••*•*•«* 


*  SUBROUTINE  BOUNDR  » 

*#*##*#»*»*****»»*#»»#*»**##»*****»»****##»*#»#*#***###*»*****#*#»*» 

SUBROUTINE  BOUNOR(C,C1  ,C2,D1 ,02,00,SG,N,NV) 

REAL  0(202,2) ,C1 ,02,01 ,D2,D0(*) ,SG(*) 

INTEGER  N,  NV 
C 

C  This  subroutine  sets  up  the  C-matrix  for  each  energy 

0  group  &  passes  it  back  to  the  main  program.  This 

0  subroutine  is  necessary  in  the  radiation  transport 

C  problem  ONLY  when  the  boundary  conditions  require  terms 

C  that  are  energy  dependent. 

C 

02  =  (0D(N-1 )*SG(N-1 ) )**0.5 
D  PRINT*," 

0  PRINT*,  'N-1  =',N-1,'DD  =\DD,'SG  =',SG 

0(1,1)  =  01 
0(2,1)  =  01 
C(NV-1 ,2)  =  02 
C(NV,2)  =  D2 

0  PRINT*,' echo  of  0-matrix  in  subroutine  BOUNDR s' 

D  DO  100  I  =  1  ,NV 

D  PRINT*,C(I ,1 ) ,0(1,2) 

D1 00  CONTINUE 
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*  SUBROUTINE  BOATA  * 

*»*«******#****»#**#*»»**#»»*»*###»****#*«*»****##*#*»»»****»*#*****#**#** 

SUBROUTINE  BOATA  (MAXG.NG.NV.CI ,C2,D1 ,D2,E1 ,E2,C) 

C 

C  this  subroutine  reads  in  boundary  data  and 

C  computes  the  boundary  matrix. 

C 

REAL  El (*),E2(*),C1 ,C2,D1 ,D2,C(202,2) 

INTEGER  MAXG,NG,I,NV 

D  PRINT*,'  entering  BOATA  subroutine  .  .  .  ' 

PRINT*, 'ENTER  Cl, C2, 01 ,02  (AS  SHOWN)' 

READ*, Cl ,C2,D1 ,02 

PRINT*, 'ENTER  ONE  VALUE  OF  El  FOR  EACH  ENERGY  GROUP' 

PRINT*, 'BEGIN  WITH  HIGHEST  GROUP  (#1),  PUTTING  EACH  ENTRY' 

PRINT*, 'ON  A  SEPARATE  LINE.' 

DO  80  I  =  1 ,NG 
READ*, El (I) 

80  CONTINUE 

PRINT*, 'NOW  REPEAT  FOR  E2' 

DO  90  I  =  1 ,NG 
READ*,E2(I) 

90  CONTINUE 
C 

C  now  construct  the  "C"  matrix  (boundary  cond.  matrix) 

C 

C(1,1)  =  Cl 
C(2,1)  =  D1 
C(NV-1,2)  =  C2 
C(NV,2)  =  02 

D  PRINT*, 'echo  the  C  matrix  in  BOATA  subroutine' 

0  DO  100  I  =  1,NV 
0  PRINT*, C(I,1 ),C(I,2) 

D100  CONTINUE 
END 
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*  SUBROUTINE  MOATA  * 

ft##*#******#*##*-*#*#*****#*##****##**##*#*##*#*####*********#***#*#*###**** 

SUBROUTINE  MDATA(MAXMAT,MAXG,NMAT,NG,XTR,XT,XS,XR) 

C 

C  this  subroutine  reads  in  material  data. 

C 

REAL  XTR(5,21 ) ,XT(5,21 ) ,XS(5,21 ,21 ) , XR(5,21 ) 

INTEGER  NG,MAXMAT,MAXG,NMAT,I,J,K 
0  PRINT*," 

D  PRINT*,'  entering  MOATA  subroutine  now  .  . 

0  PRINT*, ' ' 

PRINT*,' ENTER  THE  NUMBER  OF  DIFFERENT  MATERIALS  USED' 

READ*,NMAT 

PRINT*, 'ENTER  THE  TOTAL  NUMBER  OF  ENERGY  GROUPS' 

READ*,NG 

C  read  in  cross  sections  beginning  with  transport  and  total 

DO  20  I  =  1.NMAT 
DO  60  J  =  1,NG 

PRINT*, 'ENTER  MATERIAL* , I, ' ,  GROUP ' , J, '  CROSS  SECTIONS' 

PRINT*, 'sigma  tr,  sigma  t:  SEPARATE  BY  COMMAS.' 
READ*,XTR(I,J),XT(I,J) 

C  now  read  in  the  scatter  cross  sections 

PRINT*, 'NON  ENTER  THE  SCATTERING  CROSS  SECTIONS.' 

PRINT*, 'BEGIN  WITH  SCATTER  FROM  HIGHEST  ENERGY  GROUP' 

PRINT*, 'ANO  STEP  THROUGH  GROUPS  UNTIL  WITHIN  GROUP  SCATTER' 

PRINT*, 'EX.  1  — >IG,  2->IG,  . .  {IG-1)->IG,  IG->IG' 

DO  70  K  =  1,J 
REAO*,XS(I,J,K) 

70  CONTINUE 

C  compute  the  removal  cross  section  from  total  and  scatter 

XR(I,J)  =  XT(I,J)  -  XS(I,J,J) 

60  CONTINUE 

20  CONTINUE 
MAXG  =  NG 
END 


I 


**************************************************************************** 


*  SUBROUTINE  GOATA  * 

it*********#*##**####******##**##****#*#****#*****#*##***#*##*#******#**##*** 

SUBROUTINE  GOATA  (MAXN,MAXMAT,R,I*IAT,N,K) 

REAL  R(*) 

REAL  A,B,C 

INTEGER  NR,K,AA,MN,NS,BB,M,MAT(*),N,I 

PRINT*, 'ENTER  THE  NUMBER  OF  REGIONS' 

RE AD*, NR 

PRINT*, 'ENTER  THE  LEFT  BOUNDARY  LOCATION' 

RE AD*. A 

PRINT*, 'ENTER  THE  GEOMETRY  TYPE,  K=?' 

READ*,K 

AA=1 

DO  50  I  =  1 ,NR 

PRINT*, 'ENTER  THE  RIGHT  HAND  BOUNDARY  FOR  REGION  ' ,1 
READ*,B 

PRINT*, 'ENTER  THE  H  OF  MESH  SPACES  FOR  THIS  REGION' 

RE AO*, NS 

PRINT*, 'ENTER  THE  MATERIAL  NUMBER' 

READ*,MN 

C  =  (B-A)/NS 
BB  =  AA  +  NS 
DO  5  M  =  AA,BB 

IF  (M.EQ.AA)  THEN 
R(M)=A 
ELSE 

R(M)  =  A  +  C 
A  =  R(M) 

ENDIF 

5  CONTINUE 

DO  10  M  =  AA,  BB-1 
MAT(M)  =  MN 
10  CONTINUE 

AA  =  BB 
N  =  BB 
MAXN=N 

MAXMAT  =  BB-1 
50  CONTINUE 
END 
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MULT I -GROUP  DIFFUSION  EQUATION  CODE 
(WITH  FINITE  ELEMENT  METHOD) 

DATE  WRITTEN s  28  OCT  1984 


AUTHOR: 


L.  WAYNE  BRASURE 


ft#********##********#*##**#*#*#*******)**#**#******#********#**********#****#* 

REAL  R(  1 01  ),XTR(5,21  ),XR(5,21  ),XT(5,21 )  ,XS'S21 ,21  ),SG(100) 

REAL  C(202,2) ,8(2) ,D0( 100) ,S( 22,3) 

REAL  Cl ,C2,D1 ,02, El (21 ) ,E2(21 ) ,T(3) 

INTEGER  MAXN .MAXMAT , N , K , MAT ( 1 00 ) ,NMAT ,NC , I , NU , NG ,MAXG 

T  ( 0 )  =  1 

T ( 1  )  =  6.2832 

T(2)  =  12.5664 

MAXG  =  21 

D  PRINT*, 'entering  data  entry  phase  of  program  ...' 

D  PRINT*," 

CALL  GDATA  (MAXN, MAXMAT, R, MAT, N,K) 

NU  =  N  *  2 

CALL  MDATA (MAXMAT , MAXG , NMAT , NG , XTR , XT ,XS , XR ) 

CALL  BDATA(MAXG,NG,NU,C1 ,C2,01 ,D2,E1 ,E2,C) 

D  PRINT*," 

D  PRINT*, 'echo  of  C  matrix  within  main  program  .  .  .' 

D  PRINT*," 

D  DO  10  I  =  1,  NU 

D  PRINT*, C(I,1 ),C(I,2) 

DIO  CONTINUE 

DO  999  IG=1 ,NG 
DO  50  I  =  1,  N-1 

SG(I)  =  XR(MAT( I) ,IG) 

DD(I)  =  1 ./(3.*XTR(MAT(I ) ,IG) ) 

D  PRINT*, 'echo  SG(',I,')  =',SG(I) 

D  PRINT*, 'echo  DD(',I,')  =',DD(I) 

50  CONTINUE 

8(1 )  =  E1(IG) 

B(2)  =  E2(IG) 

CALL  B0UNDR(C,C1 ,C2,D1 ,D2,DD,SG,N,NU) 

CALL  0IFF(5,C,B,R,00,SG,NUIN,NG,XS,XR,MAT,IG,T) 

PRINT*," 

PRINT*," 

PRINT*,'  THE  SOLUTION  FOLLOWS  ===  GROUP ',IG 
PRINT*," 

PRINT*," 

0  PRINT*, S(NV-1 ,1  ) 

DO  60  I  =  1,NV,2 

IX  =  I  -  (1/2  -  0.5) 

PRINT*,'  at  R  =' ,R(IX) , 'em,  F  =',S(I,1) 

60  CONTINUE 

C  now  end  the  iteration  over  all  energy  groups  .  .  . 

999  CONTINUE 
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fluxes  and  currents  at  each  node,  it  sets  up  the  local  and 
global  matrices  for  each  group  (one  at  a  time),  performs  the 
Simpson's  approximation,  and  solves  the  series  of  matrix 
equations  (using  IMSL  subroutines)  in  order  to  arrive  at  the 
solution  for  each  group.  Each  group  calculation  is  performed 
separately,  beginning  with  the  highest  energy  group,  group 
one.  The  solution  is  stored  in  the  S  matrix,  printed  out, 
and  then  the  subroutine  moves  on  to  the  next  lower  group. 

List  of  Variables. 

NI  -  the  number  of  mesh  intervals, 

INSCATTER  -  1  if  doing  an  inscatter  calculation, 

D  -  diffusion  constant  in  current  interval, 

H  -  width  of  current  interval, 

SA,  SB,  SC  -  variables  used  in  Simpson's  approx., 

SOURCE  -  stores  the  source  term  used  in  Simpson's 
approx .  , 

XX  -  stores  the  positions  of  each  node, 

R,  Q  -  variables  used  in  solving  linear  equations 
in  conjunction  with  IMSL  subroutines. 


c 


c 

C  we  now  have  the  FLUENCE  due  to  scattering  and 

C  must  add  this  to  the  direct  fluence  at  each  node. 

C 

C  first,  we  must  calculate  the  direct  fluence  for  this 

C  energy  group  at  each  node. 

C 

CALL  DIRECT (N,T , XX, MAT, XT ,DPHI  .IGjIA.PIPHI .MIDPT ) 

PRINT*,' - ' 

PRINT*," 

PRINT*,'  Scattered  Fluence  for  Group', IG,'  is  =',S(NU-1,1) 
PRINT*,'  &  the  Direct  Fluence  =',DPHI(N) 

00  700  I  =  1  ,NV,2 

IX  =  I  -(1/2  -  0.5) 

D  PRINT*, 'S( ' ,1, ' ,1 )  =' ,S(1, 1  ) 

D  PRINT*, 'DPHIC, IX,')  ='  ,DPHI(IX) 

3(1,1)  =  S(I,1)  +  DPHI(IX) 

0  PRINT*, 'SC, I,', 1)  ='  ,S(I,1 ) 

C  put  the  solutions  for  this  energy  group  into  the 

C  "PHI"  array,  to  be  used  for  the  downscattered 

C  contribution  for  the  lower  energy  groups. 

PHI(IG.IX)  =  S(I,1) 

0  PRINT*,'  for  IX  =',IX,'  and  IG  =',IG 

D  PRINT*,'  PHI(IG.IX)  =',PHI(IG,IX),'  S(I,1)  =',S(I,1) 

700  CONTINUE 
C 

C  return  control 

C 

END 


im 


Appendix  F:  Generating  Poisson  and  Multinomial 
Reference  Values 


£ 


This  appendix  briefly  discusses  the  way  in  which  the 
reference  source  spectra  were  arrived  at  for  both  the  Poisson 
and  multinomial  versions  of  the  Bayes'  theorem  analysis. 

Both  methods  use  the  numbers  generated  by  the  "channel 
percentage  converting  program"  which  is  listed  in  this 
appendix.  This  program  converts  the  raw  transport  data 
contained  in  tables  III-l  through  III-6  into  a  form 
acceptable  to  the  Bayes'  theorem  analysis  codes  discussed  in 
appendix  H. 


m* 


Program  Algorithm  and  Use 

For  each  of  the  three  sources  at  each  of  the  ranges 
given,  the  program  inputs  the  number  of  counts  per  photon  in 
each  channel.  The  program  output  includes  the  total  counts 
per  photon  for  the  given  source  as  well  as  a  listing  of  the 
proportion  of  counts  in  each  channel  (the  sum  of  all  channels 
being  unity) . 

The  total  counts  per  photon  value  is  used  to  determine 
the  total  number  of  counts  for  each  source  at  any  given 
distance.  For  instance,  if  at  a  distance  of  50  meters, 
measurements  of  equal  time  (Poisson)  or  of  equal  counts 
(multinomial)  yield  10,000  counts  (assigned  arbitrarily), 
then  the  number  of  counts  at  all  other  ranges  can  be 
determined  for  all  sources.  For  source  A  at  100  meters,  the 
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total  counts  per  photon  (all  channels)  of  source  A  is  divided 
by  the  total  counts  per  photon  (all  channels)  of  source  A  at 
50  meters,  and  multiplied  by  the  number  of  counts  at  50 
meters.  This  yields  the  folowing:  9 . 12E-10/3 . 16E-9  *  10000  * 
2816.7.  So  that  2816.7  counts  will  be  obtained  (2817  for  the 
multinomial  analysis).  The  reference  sources  will  be 
constructed  from  this  type  of  data,  as  discussed  in  the  next 
two  sections. 

Poisson  Reference  Sources 

The  first  step  in  the  process  is  to  assign  a  level  of 
counts  at  the  50  meter  range.  Then,  using  the  method 
discussed  above,  the  number  of  counts  for  all  three  sources 
at  each  range  is  calculated.  Next,  using  the  proportion  of 
counts  in  each  channel,  the  number  of  counts  in  each  channel 
is  determined.  This  results  in  three  reference  sources  at 
each  range. 

Multinomial  Reference  Sources 

To  obtain  the  multinomial  reference  sources,  the  total 
number  of  counts  for  each  source  at  each  range  is  determined 
in  the  same  manner  as  above.  The  reference  sources  at  each 
range  consist  simply  of  the  proportions  in  each  channel 
(recall  the  form  of  the  multinomial  distribution). 

Program  Output 


The  output  from  this  program  has  been  surpressed,  since 


it  follows  readily  from  tables  III-l  through  III-6 


List  of  Variables 
integer : 

COUNTS  =  counts/photon  in  each  channel, 

N  =  source  number  (1-3), 

NN  =  scientific  notation  subroutine  variable, 
SUM  =  total  counts/photon  (all  channels), 
real : 

XX  =  scientific  notation  subroutine  variable. 
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10  REM 

20  REM  +++++++++++++++++++++++++++++++++++++♦++ 

30  REM  + 

35  REM  +  CHANNEL  %  CONVERTING  PROGRAM  + 

AO  REM  +  By  L.  liiayne  Brasure  + 

50  REM  +  December  1984  + 

55  REM  +  + 

70  REM 

80  OIM  C0UNTS(18) 

85  INPUT  "  Enter  the  source  number  =>";N 

90  INPUT  "  Enter  the  range  (in  meters)  =>"; RANGE 

100  PRINT  "  Enter  the  counts  per  photon  for  each  channel  indicated 

110  SUM  =  0 

120  FOR  I  =  1  TO  18 

130  PRINT  "  Channel  Number  ";I;"  =>" 

140  INPUT  COUNTS (I ) 

150  SUM  =  SUM  +  COUNTS(I) 

160  NEXT  I 

170  PRINT  "For  Source  Number  ";N{",  at  R  =  "{RANGE;"  meters:" 

175  PRINT 

180  PRINT  "  total  counts  per  photon  =";SUM 
185  PRINT 

190  PRINT  "  (channel  breakdown  follows)" 

200  PRINT 

210  PRINT  "  CHANNEL  NUMBER  PERCENT  OF  TOTAL  COUNTS" 

220  PRINT  "  _  _ " 

225  PRINT  "" 

230  NN  =  4 

235  REM  +++++  print  out  the  table  now  .  .  . 

240  FOR  I  =  1  TO  18 

250  COUNTS(I)  =  COUNTS(I)  /  SUM 

253  PRINT  " 

255  XX  =  I:  GOSUB  2100:  PRINT  "  "{ 

260  XX  =  COUNTS(I) 

270  GOSUB  2100 
275  PRINT 
280  NEXT  I 
290  PRINT 

300  INPUT  "  Care  for  another  run?  Enter  [1]  for  yes  =>";ANS 
310  IF  (ANS  =  1)  GOTO  85 
320  ENO 
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+  SCIENTIFIC  NOTATION  SUBROUTINE  + 
++++++++++++++++++++++++++++++++++++++++ 


2100  REM 
2110  REM 
2120  REM 
2190  REM 
2200  REM 
2210  REM 

2220  IF  NN  <  0  OR  NN  >  8  THEN  PRINT  "RANGE  ERR";:  RETURN 
2230  EX  =  0:  IF  XX  =  0  THEN  MTS  =  "0.":  GOTO  2250 

2240  MT  =  UAL  (  STR$  (  A8S  (XX))):  GOSUB  2260:  IF  NN  <  8  THEN  MT  =  MT  +  . 

5  *  10  a  (  -  NN):  GOSUB  2260 

2250  PRINT  MIDS  (XX  <  0)  +  1,1);  LEFTS  (MTS  +  "00000000" ,NN  +  2);" 

E";  MIDS  ("+-", (EX  <  0)  +  1,1);  RIGHTS  ("0"  +  STRS  (  ABS  (EX)), 2);:  RETURN 
2260  IF  MT  >  =  10  THEN  MT  =  MT  /  1Q:EX  =  EX  +  1 :  GOTO  2260 

2270  IF  MT  <  1  THEN  MT  =  MT  *  10:EX  =  EX  -  1 :  GOTO  2270 
2280  MTS  =  STRS  (MT):  IF  MIDS  (MT$,2,1)  <  >  THEN  MTS  =  MTS  + 

2290  RETURN 


Appendix  G:  Generating  Measured  Source  Spectra 


This  appendix  will  address  the  method  of  generating  the 
measured  source  spectra  during  the  analysis  section  of  the 
problem.  The  method  will  be  discussed  separately  for  each  of 
the  three  types  of  statistical  distributions  used. 

Poisson  Random  Number  Generator  ref.  (11)  &  (16) 

Algorithm  Description. 

The  Poisson  distribution  inputs  the  reference  sources  in 
the  form  of  number  of  counts  in  each  channel  (units  = 
counts),  unlike  the  multinomial  distribution  format.  One  of 
the  three  reference  sources  is  selected  as  the  source  which 
is  to  be  "measured,"  and  the  measured  source  array  is 
initially  assigned  as  this  particular  reference  source. 

The  Poisson  random  number  subroutine  first  generates  a 
random  number  between  between  0  and  1  (inclusive)  and  assigns 
the  value  to  the  variable  NUM.  An  integer  is  incremented, 
beginning  at  1,  and  is  stored  in  the  variable  N.  The 
probability  of  obtaining  this  number  N  is  calculated  using 
the  Poisson  distribution  with  the  number  of  reference  counts 
in  the  channel  serving  as  the  mean.  The  probability  is  next 
compared  with  the  random  number,  and  the  integer  N  is 
increased  until  the  probability  is  greater  than  the  random 
number,  NUM.  This  process  is  repeated  for  all  18  channels. 
List  of  Variables. 
integer : 


IG  =  channel  counter, 


N  =  number  of  measured  counts, 

U  =  number  of  counts, 
real : 

MPHI(18)  =  measured  source  array  (ends  up  an  integer), 
NUM  =  random  number  variable, 

P  =  probability  variable, 

T  =  probability  variable, 

Y  =  number  of  reference  counts  variable. 

Gaussian  Random  Number  Generator  ref.  (16) 

Algorithm  Description. 

The  Gaussian  random  number  generator  begins  in  the  same 
manner  as  does  the  Poisson  random  number  generator,  with  the 
MPHI  array  initially  equal  to  the  chosen  reference  source. 
Next,  for  each  channel  of  the  measured  source,  twelve  random 
numbers  are  generated,  since  the  uniform  random  number 
generator  approximates  the  Gaussian  distribution  with  this 
many  repetitions  (16).  Each  time  a  random  number  is 
generated,  it  is  added  to  the  number  of  counts  in  the  given 
channel  less  6.  After  all  the  numbers  have  been  generated, 
the  resultant  real  number  of  counts  is  rounded  to  an  integer 
number  of  counts. 

List  of  Variables. 
integer : 

IG  =  channel  counter, 

IMPHI(I8)  =  final  measured  source  array, 


real 


MPHI(18)  =  initial  measured  source  array, 

X  =  random  number  variable, 

Y  =  reference  counts  in  a  channel  minus  6. 

Multinomial  Random  Number  Generator  ref.  (16) 

Algorithm  Description. 

Unlike  the  Poisson  and  Gaussian  distributions,  the 
multinomial  distribution  subroutine  inputs  the  reference 
sources  as  fractions  of  counts  in  each  channel  (units  = 
counts/photon).  To  generate  a  measured  spectrum,  one  of  the 
reference  sources  is  selected  as  before.  The  fractions  in 
each  channel  are  converted  into  a  cumulative  distribution,  so 
that  the  last  channel  in  the  MPHI  array  will  contain  a  value 
of  1. 

The  multinomial  random  number  generator  then  assigns 
each  one  of  the  total  counts  available  to  one  of  the  channels 
in  the  measured  source  array  as  follows.  First,  a  random 
number  between  0  and  1  is  generated,  which  is  then  compared 
with  the  cumulative  probability  in  each  channel  of  the  MPHI 
array.  When  the  random  number  is  less  than  the  cumulative 
probability  in  one  of  the  channels,  the  count  is  added  to 
that  particular  channel.  This  process  continues  until  all 
the  available  counts  have  been  depleted,  leaving  a  new 
measured  source  spectrum. 

List  of  Variables. 

integer  : 


2 


I  =  total  count  incrementor, 

IG  =  channel  number  counter, 

MPHI(18)  =  measure  source  array, 

NC(18)  =  counting  array, 

SUM  =  total  number  of  counts  available, 

T  =  channel  number  counter, 
real : 

X  =  random  number  variable. 

Subroutine  Listings. 

The  random  number  generators  are  listed  with  the  main 


analysis  programs  in  appendix  H 


Appendix  H:  Analysis  Code  Listings 


This  appendix  describes  and  lists  the  three  codes  used 
to  generate  the  measured  source  spectra  and  perform  the 
Bayes'  theorem  analysis.  There  is  one  code  for  each  of  the 
statistical  distributions  used:  Poisson,  Gaussian  and 
multinomial.  The  first  code  listed  is  the  Poisson  version, 
and  will  be  listed  in  its  entirety.  To  obtain  the  Gaussian 
code,  two  subroutines  need  to  be  replaced:  the  RANDOM 
subroutine  and  the  POISSON  subroutine.  Subroutine  RANDOM  was 
discussed  in  appendix  F,  and  will  not  be  discussed  further. 
The  listing  of  the  two  new  subroutines  (RANDOM  and  GAUSSIAN) 
begins  on  page  132.  Similarly,  in  the  multinomial  code,  only 
those  subroutines  which  are  different  from  the  Poisson  code 
will  be  listed  and  described.  The  multinomial  listing  begins 
on  page  133. 

All  analysis  codes  were  written  in  BASIC  to  run  on  the 
Apple  II  series  of  computers,  and  minor  modifications  will 
translate  them  into  forms  which  can  be  run  on  almost  any 
microcomputer . 

The  Bayes1  with  Poisson  Program 
Main  Program. 

The  main  program  inputs  the  number  of  channels  (energy 
groups)  and  the  range  at  which  the  analysis  will  be 
performed.  It  then  calls  the  subroutines  to  solve  for  the 
Bayes'  posterior  distribution,  and  prints  out  the  results. 


The  calculations  are  performed  for  one  measured  source  at  a 
time  at  one  range. 

List  of  Variables. 

RPHI  -  stores  the  three  reference  source  spectra, 

MPHI  -  stores  the  measured  source  spectrum, 

NC  -  stores  the  number  of  counts  in  each  channel 
of  the  measured  spectrum, 

P  -  stores  P(X|Ai)  value  of  each  source, 

H$  -  string  array  used  in  reading  data 
from  data  files, 

RANGE  -  distance  from  the  source  at  which  the 
anlysis  is  performed, 

NG  -  the  number  of  channels  in  each  source, 

PI,  P2,  P3  -  comprise  the  posterior  distribution, 

NN,  XX  -  scientific  notation  variables, 

ANS  -  integer  value  1  continues  program 
execution . 

Subroutine  REFERENCE. 

This  subroutine  retrieves  the  reference  sources  from 
data  files  stored  on  disk.  It  feeds  the  data  into  the  RPHI 
array . 

Subroutine  MEASURED. 

This  subroutine  fills  up  the  MPHI  array  with  the  values 
from  whichever  source  has  been  selected. 

Subroutine  RANDOM . 

This  subroutine  generates  the  measured  source  spectrum 
using  the  Poisson  distribution.  Refer  to  appendix  F  for 
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details  on  this  subroutine. 

Subroutine  BAYES. 

This  subroutine  performs  the  Bayes'  theorem  analysis  as 
depicted  in  equation  (30)  of  chapter  IV. 

List  of  Variables. 

Al,  A2,  A3  -  the  assumed  prior  distribution 
of  the  reference  sources, 

DENOM  -  stores  the  denominator  in  equation  (30), 

XI,  X2,  X3  -  the  values  of  P(X|Ai)  returned  by 
subroutine  POISSON. 

Subroutine  POISSON. 

This  subroutine  calculates  the  values  of  P(x|Ai)  using 
the  Poisson  distribution,  as  discussed  in  chapter  IV. 

List  of  Variables. 

S  -  reference  source  counter, 

P  -  stores  the  final  values  of  P(X|Ai), 

IG  -  channel  counter, 

I  -  count  level  counter, 

FT  -  temporary  storage  variable, 

TERM  -  temporary  storage  variable. 

Scientific  Notation  Subroutine. 

This  "canned"  subroutine  formats  the  output  on  the  Apple 
He  computer  and  was  written  by  John  Baldwin  of  Erie, 
Pennsylvania. 

Data  Retrieval  Su b r outine . 


The  Baves'  with  Gaussian  Proaram 


Listing  begins  on  page  132. 

Main  Program. 

The  main  program  is  identical  to  the  Poisson  version, 
and  hence,  only  the  two  new  subroutines  will  be  discussed. 

Subroutine  RANDOM. 

This  subroutine  generates  the  measured  source  spectrum 
using  the  Gaussian  distribution.  Refer  to  appendix  F  for 
details  on  this  subroutine. 

Subroutine  GAUSSIAN. 

This  subroutine  calculates  the  values  of  P(X|Ai)  using 
the  Gaussian  distribution,  as  discussed  in  chapter  IV. 

The  Bayes'  with  Multinomial  Program 

Listing  begins  on  page  133. 

Main  Program. 

The  main  program  is  basically  the  same  as  the  Bayes' 
with  Poisson  program,  so  once  again,  only  the  new  subroutines 
will  be  discussed. 

Subroutine  MEASURED. 

The  modification  to  this  subroutine  is  to  calculate  the 
cumulative  probability  in  each  channel  of  the  array  MPHI,  in 
order  to  calculate  a  measured  spectrum  using  the  multinomial 
random  number  generator  discussed  in  appendix  F. 

Subroutine  RANDOM. 


This  subroutine  generates  the  measured  source  spectrum 


This  report  shows  that  gamma  ray  spectra  identification 
using  Bayes'  Probability  Theorem  can  be  used  to  extend  positive 
identification  ranges  when  compared  with  the  method  of  photopeak 
identification.  In  this  study,  Bayes'  Theorem  methodology  extended 
the  range  of  positive  identification  a  minimum  of  50  meters  in  a  low 
count  environment.  These  results  are  based  on  spectra  generated  using 
the  Finite  Element  Method  in  conjunction  with  Poisson  counting  statistics. 
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5550  REM  +++++  Now  factor  in  the  last  few  factorial  terms 

5560  TI  =  SUM  -  TI  -  1 

5565  IF  (TI  <  1)  GOTO  5600 

5570  FOR  I  =  TI  TO  1  STEP  -  1 

5580  P(S)  =  P(S)  *  I 

5590  NEXT  I 

5600  NEXT  S 

5610  XI  =  P(1):  PRINT  "XI  =  ";X1 

5620  X2  =  P(2):  PRINT  "X2  =  ";X2 

5630  X3  =  P(3):  PRINT  "X3  =  ";X3 

5700  RETURN 


5000  REM  +++++  a  factorial  fix-all  begins  here: 
5010  P(S)  =  P(S)  »  (SUM  -  TI) 

5020  TI  =  TI  +  1 
5030  GOTO  5240 


5000  REM 

5010  REM  ++++++++++++++++++++++++++++++++++++++++ 

5020  REM  +  SUBROUTINE  MULTINOMIAL  + 

5030  REM  ++++++++++++++++++++++++++++++++++++++++ 

5040  REM 

5050  REM  +++++  Bubble  sort  the  MPHI  array,  and  place 

5060  REM  +++++  the  RPHI  values  in  a  corresponding  order. 

5062  SWITCH  =  1 

5063  IF  (SWITCH  =  0)  GOTO  5080s  REM  (go  to  end  statement). 

5064  SWITCH  =  0:  REM  (0  if  sort  is  complete) 

5065  FOR  I  =  2  TO  NG 

5066  IF  (MPHI(I)  <  MPHI (I  -  1))  GOTO  5068 

5067  GOTO  5078 

5068  TEMP  =  MPHI(I) 

5069  MPHI ( I )  =  MPHI (I  -  1) 

5070  MPHI (I  -  1)  =  TEMP 

5071  REM  +++++  Now,  the  fractions  must  move  WITH  the  MPHI  sort 

5072  FOR  IV  =  1  TO  3 

5073  T9  =  RPHI (IV, I) 

5074  RPHI (IV, I)  =  RPHI (IV, I  -  l) 

5075  RPHI (IV, I  -  1)  =  T9 

5076  NEXT  IV 

5077  SWITCH  =  1 

5078  NEXT  I 

5079  GOTO  5063 


+++++  Bubble  sort  is  now  complete  .  .  . 

+++++  Now,  calculate  the  P(X(source  1),  etc. 
S  =  1  TO  3 
=  1 

MPHI(NG) 

1 


5080  REM 
5090  REM 
5100  FOR 
5150  P(S) 

5200  IM  = 

5210  JL  =  1 sTI 
5220  P(S)  =  SUM 
5230  FOR  I  =  1  TO  IM 

5240  IF  (TI  <  SUM  AND  P(S)  <  1E30)  GOTO  5800 

5250  IF  ((MPHI(JL)  -  I)  <  0)  GOTO  5300 

5260  GOTO  5500 

5300  JL  =  JL  +  1 

5400  GOTO  5250 

5500  FOR  J  =  JL  TO  NG 

5510  P(S)  =  P(S)  »  RPHI(5,J)  /  I 

5520  NEXT  J 

5530  NEXT  I 


2015  REM  +  SUBROUTINE  MEASURED  + 

2020  REM  ++++++++++++++++++++++++++++++++++++++++ 

2025  REM 

2030  PRINT  s  PRINT  "MEASURED  SOURCE  SPECTRUM":  PRINT  " - 

_ t! 

2040  PRINT  "  at  RANGE;"  cms.":  PRINT 

2050  PRINT  "enter  the  applicable  ref.  source  (1,2,or3)  PRIN 
2055  PRINT  "  (the  cumulative  prob.  will  be  calculated)" 

2057  INPUT  I 

2060  FOR  IG  =  1  TO  NG 

2065  FRACT(IG)  =  PHI(I.IG) 

2081  MPHI(IG)  =  FRACT(IG) 

2085  IF  (IG  >  1)  THEN  MPHI(IG)  =  MPHI(IG)  +  MPHI(IG  -  1) 

2090  PRINT  "  MPHI(";IG;")  =";MPHI(IG) FRACT(IG)  =  ";FRACT(IG) 
2100  NEXT  IG 
2110  RETURN 
3000  REM 

3020  REM  +  SUBROUTINE  RANDOM  + 

3030  REM  +++++++++++++++++++++++++++++++++++++ 

3040  REM 

3041  REM  +++++  initialize  the  counter  array: 

3042  FOR  T  =  1  TO  18:NC(T)  =  0:  NEXT  T 
3060  FOR  1=1  TO  SUM 

3065  X  =  RND  (1) 

3070  FOR  IG  =  1  TO  NG 
3080  Y  =  X  -  MPHI(IG) 

3090  IF  (Y  <  0)  THEN  GOTO  3100 
3095  NEXT  IG 
3100  NC(IG)  =  NC(IG)  +  1 
3110  NEXT  I 

3120  PRINT  :  PRINT  "  Measured  source  spectrum  follows:":  PRINT 

3125  FOR  IB  =  1  TO  3:  PRINT  CHR$  (7):  NEXT  IB 

3130  FOR  IG  =  1  TO  NG 

3140  PRINT  "  channel  §  ";IG;"  =  ”;NC(IG) 

3150  MPHI(IG)  =  NC(IG) 

3155  NEXT  IG 
3200  RETURN 


3000  REM 

3010  REM  +++++++++++++++++++++++++++++++++++++ 

3020  REM  +  SUBROUTINE  RANDOM  + 

3030  REM  +++++++++++++++++++++++++++++++++++++ 

3040  REM 

3050  REM  +++++  This  subroutine  generates  random  numbers 

3060  REM  +++++  FOLLOWING  A  GAUSSIAN  DISTRIBUTION  .  .  . 

3070  FOR  IG  =  1  TO  NG 

3000  Y  =  MPHI(IG)  -  6.0 

3090  FOR  J  =  1  TO  12 

3100  X  =  RND  (1) 

3200  Y  =  Y  +  X 
3210  NEXT  J 
3220  MPHI(IG)  =  Y 

3230  IF  (MPHI(IG)  <  0)  THEN  MPHI(IG)  =  0 
3240  NEXT  IG 

3300  REM  +++++  Echo  the  measured  source  spectrum:” 

3310  PRINT  CHR$  (7):  PRINT  CHR$  (7) 

3320  PRINT  "Measured  source  spectrum  follows:":  PRINT 

3330  FOR  IG  =  1  TO  NG 

3335  IMPHI(IG)  =  INT  (MPHI(IG)  +  0.5) 

3340  PRINT  "channel  ";IG;"  =  ";IMPHI(IG) 

3350  NEXT  IG 
3360  PRINT 
3400  RETURN 


5000  REM 

5010  REM  ++++++++++++++++++++++++++++++++++++++++ 

5020  REM  +  SUBROUTINE  GAUSSIAN  + 

5030  REM  ++++++++++++++++++++++++++++++++++++++++ 

5040  REM 

5100  FOR  S  =  1  TO  3 

5110  P(S)  =  1 

5120  DUMMY  =  1 

5130  FOR  IG  =  1  TO  NG 

5140  IF  (IMPHI(IG)  >  0)  GOTO  5200 

5150  ZT  =  1 

5160  GOTO  5300 

5200  ZT  =  IMPHI(IG) 

5300  DUMMY  =  EXP  (  -  0.5  *  ((IMPHI(IG)  -  RPHI(S.IG))  /  ZT  1  0.5)  ’  2)  /  ( 
2  *  3.142593  *  ZT)  J  0.5 
5310  P(S)  =  P(S)  *  OUMMY 
5320  NEXT  IG 
5330  NEXT  S 

5440  XI  =  P(1 ) :  PRINT  "XI  =  ";X1 

5450  X2  =  P(2) :  PRINT  "X2  =  ";X2 

5460  X3  =  P(3):  PRINT  "X3  =  ";X3 

5500  RETURN 


132 


1  TO  3 


REM 

REM  +  SUBROUTINE  POISSON  + 

REM  +++++++++++++++++++++++++++++++++-*•++++++ 

REM 

FOR  S  =  1  TO  3 
P(S)  =  1 
FT  =  1 
DUMMY  =  1 
FOR  IG  =  1  TO  NG 

REM  +++++  calculate  factorial  in  numerator 
FT  =  1 

FOR  I  =  1  TO  MPHI(IG) 

FT  =  FT  *  I 
NEXT  I 

IF  (RPHI(S.IG)  >  0)  GOTO  5300 
TERM  =  1 
GOTO  5400 

TERM  =  RPHI(S.IG)  1  MPHI(IG) 

DUMMY  =  EXP  (  -  RPHI(S,IG) )  »  TERM  /  FT 
P(S)  =  P(S)  *  OUMMY 
NEXT  IG 
NEXT  S 

XI  =  P(1 ):  PRINT  "XI  =  ";X1 

X2  =  P(2):  PRINT  "X2  =  ";X2 

X3  =  P(3):  PRINT  "X3  =  ";X3 

RETURN 
REM 

REM  ++++++++++++++++++++++++++++++++++++++++ 

REM  +  SCIENTIFIC  NOTATION  SUBROUTINE  + 


IF  NN  <  0  OR  NN  >  8  THEN  PRINT  "RANGE  ERR";:  RETURN 
EX  =  0:  IF  XX  =  0  THEN  MT$  =  "0.":  GOTO  6150 

MT  =  UAL  (  STR$  (  ABS  (XX))):  GOSUB  6160:  IF  NN  <  8  THEN  MT  =  MT  +  . 

5  *  10  1  (  -  NN):  GOSUB  6160 

PRINT  MID$  ("+-", (XX  <  0)  +  1,1);  LEFTS  (MT$  +  "00000000", NN  +  2);" 

E";  MIDS  ("+-", (EX  <  0)  +  1,1);  RIGHTS  ("0"  +  STR$  (  ABS  (EX)), 2);:  RETURN 
IF  MT  >  =  10  THEN  MT  =  MT  /  10:EX  =  EX  +  1 :  GOTO  6160 

IF  MT  <  1  THEN  MT  =  MT  *  10:EX  =  EX  -  1 :  GOTO  6170 
MTS  =  STRS  (MT):  IF  MIDS  (MTS, 2,1)  <  >  "."  THEN  MTS  =  MT$  + 

RETURN 

REM 

REM  *  DATA  RETRIEVAL  SUBROUTINE  + 

REM  +♦++++++++++++++++++++++++++++++++++++++++ 

GET  C$ 

IF  C$  =  R$  THEN  RETURN 
H$( IG)  =  H$(IG)  +  C$ 

GOTO  7040 


3000  REM 

3010  REM  +++++++++++++++++++++++++++++++++++++ 

3020  REM  +  SUBROUTINE  RANDOM  + 

3040  REM 

3050  REM  +++++  This  subroutine  generates  random  numbers 

3060  REM  +++++  following  a  Poisson  distribution  .  .  . 

3070  FOR  IG  =  1  TO  NG 
3080  Y  =  MPHI(IG) 

3100  IF  (Y  <  =  0)  GOTO  3300 

3110  NUM  =  RNO  (1) 

3120  N  =  0 

3130  T  =  EXP  (  -  Y) 

3140  P  =  T 

3150  REM  +++++  The  00  UNTIL  loop  begins  here: 

3160  IF  (NUM  <  =  P)  GOTO  3400 

3170  N  =  N  +  1 

3190  U  =  N 

3200  T  =  T  *  Y  /  U 

3210  P  =  P  +  T 

3220  GOTO  31 60 

3300  REM  +++++  else,  N  =  Os 

3310  N  =  0 

3400  REM  +++++  assign  new  value  of  MPHI(IG) 

3410  MPHI(IG)  =  N 
3500  NEXT  IG 

3510  PRINT  CHR$  (7)s  PRINT  CHR$  (7) 

3600  PRINT  s  PRINT  "  the  measured  spectrum  f ollous : " :  PRINT 

3610  FOR  IG  =  1  TO  NG 

3620  PRINT  "  channel  ";IG;"  =  ";MPHI(IG) 

3630  NEXT  IG 
3700  RETURN 
4000  REM 

4020  REM  +  SUBROUTINE  BAYES  + 

4040  REM 

4050  REM  +++++  First,  assign  the  prior  Distribution  values 

4060  REM  +++++  of  the  three  reference  sources  .  .  , 

4070  A1  =  0.3333sA2  =  0.3333:A3  =  0.3333 

4080  REM  +++++  Call  the  statistics  subroutine  to 

4090  REM  +++++  determine  the  values  of  P(X|source  1), 

4100  REM  +++++  P  (Xjsource  2),  and  P  (X (source  3). 

4110  GOSUB  5000 

4120  REM  +++++  Calculate  the  denominator  term  of  Bayes'  Thm 
4130  OENOM  =  XI  *  A1  +  X2  »  A2  +  X3  *  A3 
4135  IF  (DENOM  =  0)  THEN  PRINT  "  OENOM  =  Os":  GOTO  4200 
4140  REM  +++++  calculate  the  likelihoods  at  the  reference 
4150  REM  +++++  range  from  the  source. 

4155  FOR  IB  =  1  TO  3s  PRINT  CHR$  (7):  NEXT  IB 

4160  PI  =  (XI  *  A1)  /  OENOM 

4165  P2  =  (X2  »  A2)  /  DENOM 

4170  P3  =  (X3  *  A3)  /  DENOM 

4200  RETURN 


1010  REPO  +  SUBROUTINE  REFERENCE  + 

1030  REN 

1040  PRINT  :  PRINT  "REFERENCE  SOURCES":  PRINT  " - 

1050  PRINT  "  at  a  range  of  "jRANGE;"  cms" 

1055  PRINT 

1056  PRINT  "  now  acquiring  necessary  data  ..." 

1057  0$  =  CHR$  (4):  REM  +++++  (ESCAPE  CHARACTER) 

1058  R$  =  CHR$  (13):  REM  +++++  (RETURN  CHARACTER) 

1059  PRINT  D$;"OPEN  SOURCE  DATA" 

1060  PRINT  D$;"READ  SOURCE  DATA" 

1070  FOR  I  =  1  TO  3 

1080  REM  +++++  loop  through  all  three  sources 

1090  FOR  IG  =  1  TO  NG 

1100  REM  +++++  loop  through  all  the  channels 

1110  GOSUB  7000 

1120  RPHI(I.IG)  =  UAL  (H$(IG)) 

1130  H$( IG)  =  "" 

1140  NEXT  IG 
1150  NEXT  I 
1155  PRINT  R$ 

1160  PRINT  D$; "CLOSE  SOURCE  DATA" 

1170  PRINT  "ECHO  DATA  ..." 

1180  PRINT  :  FOR  I  =  1  TO  3 
1185  PRINT  "  SOURCE  §  ";I 
1190  FOR  IG  =  1  TO  NG 

1200  PRINT  "  channel  ";IG;"  =  ";RPHI(I,IG) 

1210  NEXT  IG 
1220  NEXT  I 
1500  RETURN 
2000  REM 

2015  REM  +  SUBROUTINE  MEASURED  + 

2020  REM  ++++++++++++++++++++++++++++++++++++++++ 
2025  REM 

2030  PRINT  :  PRINT  "MEASURED  SOURCE  SPECTRUM":  PRINT  " 

_ It 

2040  PRINT  "  at  "jRANGEj"  cms.":  PRINT 

2050  PRINT  "enter  the  applicable  ref.  source  (1,2,or3) 

2055  PRINT  "  (the  cumulative  prob.  will  be  calculated) 

2057  INPUT  I 

2060  FOR  IG  =  1  TO  NG 

2065  FRACT(IG)  =  RPHI(I,IG) 

2081  MPHI(IG)  =  FRACT(IG) 

2090  PRINT  "  MPHI(";IG;")  =";MPHI(IG) 

2100  NEXT  IG 
2110  RETURN 


30  REM  +  + 

40  REM  +  BAYES  PROBABILITY  ANALYSIS  + 

50  REM  +  (using  the  Poisson  distribution)  + 

60  REM  +  By:  L.  Wayne  Brasure  + 

70  REM  +  December  1984  + 

80  REM  +  + 

90  REM  ++++++++++++++++++++++++++++++♦+++++++++ 

100  REM 

110  DIM  RPHI(3,18):  OIM  MPHI(18):  DIM  NC(18):  DIM  FRACT(IB) 

111  DIM  P(3):  DIM  H$(18) 

115  PRINT  "random  number  seed  =  RNO  (  -  7) 

118  PRINT  :  PRINT  :  HOME 

120  PRINT  "BAYES  PROBABILITY  ANALYSIS" 

130  PRINT  "(the  multinomial  version  w/no  background  .  .  .)" 

140  PRINT  :  PRINT 

150  INPUT  "Enter  the  total  number  of  channels  =>  ";NG 
160  INPUT  "Enter  the  distance  from  the  source  (cms)  =>  "; RANGE 
170  REM  +++++  enter  the  reference  source  spectra  at  this  range 

180  GOSUB  1000 

190  REM  +++++  enter  the  measured  source  spectrum  9  this  range 
200  GOSUB  2000 

210  REM  +++++  Calculate  measured  source  spectrum 
215  REM  +++++  using  POISSON  random  number  generator 
220  GOSUB  3000 

230  REM  +++++  Calculate  the  Bayesian  Posterior  Distribution 

240  GOSUB  4000 

245  IF  (DENOM  =  0)  GOTO  360 

250  REM  +++++  Print  out  the  Posterior  Distribution 
260  NN  =  4 
265  PRINT 

270  PRINT  "  SOURCE  NUMBER  PROB  (SOURCE(i)iX)" 

280  PRINT  "  _  _  " 

290  PRINT 

300  PRINT  "A  "{ 

310  XX  =  PI:  GOSUB  6000:  PRINT 

320  PRINT  "  B 

330  XX  =  P2:  GOSUB  6000:  PRINT 

340  PRINT  ”  C  "; 

350  XX  =  P3:  GOSUB  6000:  PRINT 

360  PRINT  :  INPUT  "  Another  run  perhaps?  i  1=yes  J  =>  ";ANS 

370  IF  (ANS  <  >1)  GOTO  995 

371  INPUT  "...  with  the  same  setup?"; ANS 

372  IF  (ANS  =  1)  GOTO  500 

380  INPUT  "...  with  just  a  new  MEASURED  source?  =>  ";ANS 
390  IF  (ANS  =  1)  GOTO  190 

400  PRINT  "  okay,  from  the  beginning  then  .  .  ." 

410  GOTO  140 

500  FOR  IG  =  1  TO  NG 

505  MPHI(IG)  =  FRACT(IG) 

520  NEXT  IG 
530  GOTO  210 

995  PRINT  :  PRINT  "Very  well  then,  good  day  sir." 

999  END 


using  the  multinomial  distribution, 
details  on  this  subroutine. 


Refer  to  appendix  F  for 


Subroutine  MliLTINOMI AL . 

This  subroutine  calculates  the  values  of 
the  multinomial  distribution,  as  discussed  in 
List  of  Variables. 

SWITCH  -  on/off  determining  variable  for 
bubble  sort, 

T9  -  a  temporary  storage  variable  for  the 
bubble  sorting  process, 

JL,  TI  -  incremental  markers  for  the 
multinomial  calculations, 

SUM  -number  of  the  total  fixed  number  of  counts 
available . 


P(X|Ai)  using 
chapter  IV. 

the 


END 
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